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Abstract 


ODRPACK  is  a software  package  for  weighted  orthogonal  distance  regression,  i.e.,  for 
finding  the  parameters  that  minimize  the  sum  of  the  squared  weighted  orthogonal  dis- 
tances from  a set  of  observations  to  the  curve  or  surface  determined  by  the  parameters. 
It  can  also  be  used  to  solve  the  nonlinear  ordinary  least  squares  problem.  The  proce- 
dure has  application  to  curve  and  surface  fitting,  and  to  measurement  error  models  in 
statistics.  ODRPACK  can  handle  both  explicit  and  implicit  models,  and  wiU  easily  ac- 
commodate complex  and  other  types  of  multiresponse  data.  The  algorithm  implemented 
is  an  efficient  and  stable  trust  region  Levenberg-Marquaxdt  procedure  that  exploits  the 
structure  of  the  problem  so  that  the  computational  cost  per  iteration  is  equal  to  that 
for  the  same  type  of  algorithm  applied  to  the  nonlinear  ordinary  least  squares  problem. 
The  package  allows  a general  weighting  scheme,  provides  for  finite  difference  derivatives, 
and  contains  extensive  error  checking  and  report  generating  facilities. 


Keywords:  orthogonal  distance  regression;  measurement  error  models;  errors  in  vari- 
ables; nonlinear  least  squares. 


Categories:  G2E,  IlBl 
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2.01  (06-19-92)  ODRPACK  2.01  corrects  two  minor  errors.  The  first  affected  the  deriva- 
tive checking  procedure,  and  the  second  affected  the  values  specified  by 
NDIGIT.  The  solutions  found  by  ODRPACK  2.00  were  not  affected  in 
any  way  by  the  first  of  these  errors,  and  were  not  seriously  affected  by 
the  second. 

Version  2.01  also  encorporates  several  other  minor  modifications  that 
further  improve  the  code.  The  most  significant  of  these  modifications 
changes  how  the  scale  values  are  used  within  the  solution  procedure. 
This  change  should  increase  the  robustness  of  the  procedure  for  poorly 
scaled  problems,  and  will  result  in  slightly  different  computed  results 
for  all  problems.  Another  of  these  changes  is  to  the  report  generated  by 
the  derivative  checking  procedure.  This  report  now  includes  the  user 
supplied  derivative  value,  as  weU  as  the  relative  difference  between  the 
user  supplied  value  and  the  finite  difference  value  it  was  checked  against. 
This  should  help  users  assess  the  validity  of  the  checking  procedure 
results. 

The  “ODRPACK  2.01  User’s  Reference  Guide”  is  essentially  the  same 
as  the  version  2.00  Guide  (NISTIR  89-4103,  Revised).  The  major  dif- 
ferences are  that 

• the  minimum  length  of  array  WORK  has  been  increased  by  (p  -|-  m)q 
locations,  and 

• the  sample  programs  and  output  shown  in  Chapter  2,  and  the  discus- 
sion of  the  derivative  checking  results  in  Chapter  5 have  been  modified 
to  reflect  the  above  mentioned  changes. 

2.00  (03-04-92)  ODRPACK  2.00  adds  several  new  features  to  those  available  in  version 
1.8  and  earlier  versions. 

• With  version  2.00,  ODRPACK  can  now  easily  accommodate  complex 
and  other  types  of  multiresponse  data,  i.e.,  data  where  each  observe^ 
tion  is  multidimensional. 
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• It  can  handle  implicit  as  well  as  explicit  models. 

• The  weighting  scheme  has  been  enhanced  to  allow  for  instances  when 
the  components  of  a multidimensional  observation  are  correlated. 

• A fcLcility  for  computing  central  finite  difference  derivatives  has  been 
added. 

• The  amount  of  work  space  required  for  the  computations  has  been 
reduced  for  most  problems. 

Because  of  these  new  features,  the  argument  lists  for  the  ODRPACK 
2.00  user  callable  subroutines,  including  the  user  supplied  subroutine 
FCN  that  now  comprises  the  functionality  of  subroutines  FUN  and  JkC, 
are  not  the  same  as  for  earlier  versions.  The  interpretations  of  some 
arguments  have  changed,  and  new  arguments  have  been  added.  In 
addition,  the  arguments  have  been  reordered,  with  arguments  having 
similar  purpose  now  grouped  together.  Those  who  have  been  using  an 
earlier  release  of  ODRPACK  should  be  especially  careful  to  examine 
the  new  calling  sequences,  and  the  definitions  for  subroutine  arguments 
FCN,  JOB,  WE  and  WD. 

1.80  (12-18-90)  ODRPACK  1.80,  a test  release  with  only  limited  distribution,  differs 
from  prior  releases  in  several  respects. 

• The  number  of  function  evaluations  required  to  find  the  solution  has 
been  reduced,  and  the  “restart”  facility  has  been  modified  to  better 
accommodate  cases  where  the  user  supplied  subroutines  FUN  and/or 
JAC  are  very  time  consuming. 

• The  printed  reports  have  been  redesigned  to  identify  parameters  that 
induce  rank  deficiency,  and  to  include  a 95%  confidence  interval  for 
the  estimated  parameters. 

• Several  enhancements  have  been  added  to  the  covariance  matrix  com- 
putations, including  the  option  of  constructing  the  covariance  matrix 
without  incurring  any  additional  derivative  evaluations. 

• The  finite  difference  approximation  to  the  Jacobian  matrix  has  been 
improved. 

1.71  (07-27-89)  ODRPACK  1.71  corrects  an  error  in  the  code  that  performs  the  compu- 
tation of  finite  difference  derivatives  with  respect  to  the  errors  A when 
m > 2 and  the  default  value  of  IFIXX  is  invoked.  (The  default  value  of 
IFIXX  is  invoked  when  IFIXX  (1,1)  is  set  to  a negative  value  or  when 
ODRPACK  routines  DODR  or  SODR  are  called.)  This  error  could  result 
in  incorrect  “fixing”  of  the  explanatory  variables,  which  woixld  affect 
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the  final  solution.  Such  “fixing”  could  be  detected  by  observing  the 
presence  of  estimated  values  for  A that  are  identically  zero.  The  error 
could  go  undetected  by  the  user,  however,  if  the  values  of  A were  not 
examined  after  the  fit. 
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Preface 


ODRPACK  is  a portable  collection  of  ANSI  Fortran  77  subroutines  for  fitting  a model 
to  data.  It  is  designed  primarily  for  instances  when  all  of  the  variables  have  significant 
errors  (see,  e.g.,  [Fuller,  1987]  and  [Boggs  and  Rogers,  1990]),  implementing  a highly  ef- 
ficient algorithm  for  solving  the  weighted  orthogonal  distance  regression  problem  [Boggs 
et  ai,  1987  and  1989],  i.e.,  for  noinimizing  the  sum  of  the  squares  of  the  weighted  orthog- 
onal distajices  between  each  data  point  and  the  curve  described  by  the  model  equation. 
It  can  also  be  used  to  solve  the  ordinary  least  squares  problem  where  all  of  the  errors 
are  attributed  to  the  observations  of  the  dependent  variable. 

ODRPACK  is  designed  to  handle  many  levels  of  user  sophistication  and  problem  diffi- 
culty. 

• It  is  easy  to  use,  providing  two  levels  of  user  control  of  the  computations,  ex- 
tensive error  handling  facilities,  optional  printed  reports,  and  no  size  restrictions 
other  than  effective  machine  size. 

• It  can  handle  implicit  as  well  as  explicit  models,  and  can  accommodate  com- 
plex and  other  types  of  multiresponse  data,  i.e.,  data  where  each  observation  is 
multidimensional. 

• The  necessary  derivatives  (Jacobian  matrices)  are  approximated  numerically  if 
they  are  not  supplied  by  the  user,  and  the  correctness  of  user  supplied  derivatives 
can  be  verified  by  the  derivative  checking  procedure  provided. 

• Both  weighted  and  unweighted  analysis  can  be  performed. 

• Subsets  of  the  unknowns  can  be  treated  as  constants  with  their  values  held  fixed 
at  their  input  values,  allowing  the  user  to  examine  the  results  obtained  by  esti- 
mating subsets  of  the  unknowns  of  a general  model  without  rewriting  the  model 
subroutine. 

• The  covariance  matrix  and  the  standard  errors  of  the  model  parameter  estimators 
are  optionally  provided. 

• The  ODRPACK  scaling  algorithm  automatically  compensates  for  poorly  scaled 
problems,  i.e.,  problems  with  model  parameters,  and/or  unknown  errors  in  the 
explanatory  variables  that  vary  widely  in  magnitude. 
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• The  trust  region  Levenberg-Marquardt  algorithm  implemented  by  ODRPACK 
[Boggs  et  al.,  1987  and  1989]  has  a computational  effort  per  step  that  is  of  the 
same  order  as  that  required  for  ordinary  least  squares,  even  though  the  number  of 
unknowns  estimated  in  the  orthogonal  distance  regression  problem  is  the  number 
of  unknown  model  parameters  plus  the  number  of  explanatory  variables,  while  the 
number  of  unknowns  estimated  in  the  ordinary  least  squares  problem  is  simply 
the  number  of  unknown  model  parameters. 

• The  code  is  portable  and  is  easily  used  with  other  Fortran  subroutine  libraries. 

Computer  facilities  for  the  ODRPACK  project  have  been  provided  by  the  National  In- 
stitute of  Standards  and  Technology  (NIST),  Applied  and  Computational  Mathematics 
Division,  and  we  gratefully  acknowledge  their  support.  Machine  dependent  constants 
are  supplied  using  subroutines  based  on  the  Bell  Laboratories  “Framework  for  a Portable 
Library”  [Fox  et  a/.,  1978a].  We  have  also  used  subroutines  from  LINPACK  [Dongarra 
et  ai,  1979],  and  from  the  “Basic  Linear  Algebra  Subprograms  for  Fortran  Usage” 
[Lawson  et  ai,  1979].  The  code  that  computes  the  t-values  used  in  constructing  the 
confidence  intervals  is  from  DATAPAC  [Filliben,  1977],  and  the  code  that  checks  user 
supplied  derivatives  was  adapted  from  STARPAC  [Donaldson  and  Tryon,  1986]  using 
algorithms  described  in  [Schnabel,  1982]. 

We  appreciate  the  comments  and  suggestions  we  have  received  regarding  earlier  versions 
of  ODRPACK  and  its  documentation.  In  particiilar,  many  of  the  improvements  in  ver- 
sion 2.01  are  the  result  of  the  authors’  collaboration  with  Paul  D.  Domich,  NIST  Applied 
and  Computational  Mathematics  Division,  and  with  David  A.  Vorp,  University  of  Pitts- 
burgh Department  of  Surgery.  We  wish  to  especially  thank  Paul  D.  Domich  for  many 
productive  discussions.  While  it  is  not  possible  to  list  everyone  else  who  has  contributed 
to  ODRPACK,  we  would  like  to  thank  Vincent  D.  Arp  (formerly  of  the  NIST  Chemical 
Engineering  Science  Division),  Hariharan  K.  Iyer  (Colorado  State  University),  Susan- 
nah B.  Schiller  (NIST  Statistical  Engineering  Division),  Bernard  Thiesse  (Ecole  Na- 
tionale  Superieure  d’Electrtechnique,  d’Electronique  d’Informatique  et  d’Hydraulique, 
Toulouse,  France),  and  Eric  J.  Vanzura  (NIST  Electromagnetic  Fields  Division),  for 
providing  us  with  data  sets  that  were  invaluable  in  testing  the  new  release. 


Paul  T.  Boggs 
Richard  H.  Byrd 
Janet  E.  Rogers 
Robert  B.  Schnabel 


June  1992 
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1.  Getting  Started 


ODRPACK  is  a portable  collection  of  ANSI  Fortran  77  subroutines  for  fitting  a model 
to  data.  It  is  designed  primarily  for  instances  when  all  of  the  variables  have  significant 
errors  (see,  e.g.,  [Fuller,  1987]  and  [Boggs  and  Rogers,  1990]),  implementing  a highly  ef- 
ficient algorithm  for  solving  the  weighted  orthogonal  distance  regression  problem  [Boggs 
et  al.f  1987  and  1989].  It  can  also  be  used  to  solve  the  ordinary  least  squares  problem, 
however,  where  all  of  the  errors  are  attributed  to  the  observations  of  only  one  of  the 
variables. 

We  suggest  that  first  time  users  of  ODRPACK  begin  by  at  least  skimming  this  chapter, 
and  then  proceed  to  the  sample  programs  given  in  §2.C,  modifying  them  as  necessary 
to  solve  their  particular  problem.  The  sample  programs  are  distributed  in  machine 
readable  form  with  the  ODRPACK  release,  and  thus  should  be  available  for  use  as 
templates  without  requiring  that  the  user  enter  the  code.  The  subroutine  arguments 
used  in  these  examples  are  defined  in  detail  in  §2.B.ii. 

Chapter  1 presents  information  that  is  especially  important  to  first  time  users  of  ODR- 
PACK. Users  are  directed  to  §1.A  for  a brief  description  of  the  orthogonal  distance 
regression  problem.  This  section  introduces  notation  and  provides  background  material 
for  understanding  the  remainder  of  the  documentation.  The  algorithm  itself  is  described 
briefly  in  §1.B,  and  in  greater  detail  in  [Boggs  et  ai,  1987  and  1989].  Options  avail- 
able for  specifying  the  problem  are  briefly  discussed  in  §1.C,  and  ODRPACK’s  report 
generation  facility  is  discussed  in  §1.D.  The  need  for  good  starting  values  for  the  un- 
known parameters  of  the  model  is  explained  in  §1.E,  and  the  use  of  weights  in  §1.F. 
Finally,  §1.G  describes  two  features  of  ODRPACK  that  simplify  the  user  interface  with 
the  package. 

Chapter  2 describes  the  use  of  ODRPACK  in  detail.  Most  experienced  users  of  ODR- 
PACK will  only  need  the  information  in  this  chapter.  The  declaration  and  call  state- 
ments for  each  user  callable  ODRPACK  subroutine  are  given  in  §2.A,  the  subroutine 
cirguments  are  defined  in  §2.B,  and  sample  programs  are  shown  in  §2.C. 

The  last  three  chapters  provide  auxiliary  information.  Chapter  3 is  generally  only  needed 
by  users  with  very  time  consuming  functions;  it  explains  how  ODRPACK’s  features  can 
be  exploited  in  order  to  minimize  cost,  and  to  reduce  the  possibility  of  losing  results 
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because  limits  imposed  by  their  computer  system  were  reached.  Chapter  4 describes 
the  details  of  certain  ODRPACK  computations,  and  Chapter  5 describes  how  the  user 
can  extract  computed  results  from  the  work  vectors. 

l.A.  Notation  and  Problem  Definition 

The  method  of  least  squares  is  often  used  to  find  the  parameters  G 3?**  of  a mathe- 
matical model  that  defines  a relationship  between  variables  that  are  subject  to  errors. 
When  errors  occur  in  more  them  one  of  the  variables,  then  the  relationship  between 
the  variables  is  frequently  referred  to  as  a measurement  error  model  or  an  errors  in 
variables  problem.  The  computational  problem  associated  with  finding  the  maximum 
likelihood  estimators  of  the  parameters  of  such  models  is  known  as  orthogonal  distance 
regression.  These  models  are  discussed  in  [Fuller,  1987],  and  [Seber  and  Wild,  1989]. 
The  orthogonal  distance  regression  procedure  implemented  in  ODRPACK  is  reviewed 
in  [Boggs  and  Rogers,  1990].  That  publication  also  summarizes  the  results  of  several 
simulation  studies  that  compare  orthogonal  distcmce  regression  results  to  those  obtained 
using  ordinary  least  squares,  where  the  errors  are  assumed  to  occur  in  only  one  of  the 
variables.  (See  also  [Boggs  et  al,  1988].) 

The  model  / that  defines  the  relationship  between  the  variables  can  be  either  linear  or 
nonlinear  in  its  parameters  p.  Sometimes  one  of  the  observed  variables  is  distinguished 
as  being  a response  that  is  dependent  upon  the  remaining  variables,  which  are  commonly 
called  the  explanatory,  regressor  ot  independent  variables.  In  these  cases,  the  explana- 
tory variables  are  often  used  to  predict  the  behavior  of  the  response  variable.  In  other 
cases,  however,  there  is  no  such  distinction  among  the  variables. 

We  say  that  there  is  an  explicit  relationship  / between  the  variables  {x,y)  if 

y«/(z;y9),  (1.1) 

where  y denotes  the  response  variable,  x denotes  the  explanatory  variables,  and  y is 
assumed  to  be  only  approximately  equal  to  f{x\P)  because  of  the  measurement  errors 
in  y and  possibly  x.  When  both  x and  y are  observed  with  error,  then  the  parameters 
P can  be  found  by  orthogonal  distance  regression.  If  only  y is  subject  to  measurement 
error,  and  x is  observed  without  error,  then  the  parameters  of  such  an  explicit  model 
can  be  obtained  using  ordinary  least  squares  procedures. 

Example:  The  model  shown  in  §2.C.i, 

Vi  « fi{xi\P)  = /5i  + /32[e^3*‘ - 1]^ 

for  i = l,...,n,  which  is  example  3.2.2  on  pages  230-238  of  [Fuller,  1987],  is 
an  explicit  orthogonal  distance  regression  model.  In  this  case,  x is  the  percent 
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saturation  of  nitrogen  gas  in  a brine  solution  forced  into  the  pores  of  sandstone, 
and  y is  the  observed  compressional  wave  velocity  of  ultrasonic  signals  propagated 
through  the  sandstone.  Here,  as  in  most  explicit  models,  the  explanatory  variable 
X is  used  primarily  to  predict  the  behavior  of  the  response  variable  y.  □ 

If  the  relationship  between  the  variables  is  expressed  implicitly,  then 

/(x;/3)  w 0 , 

where  here  we  denote  all  of  the  variables  by  x.  Although  one  of  these  variables  may 
still  be  dependent  on  the  others,  when  the  model  is  implicit  we  do  not  assume  that  such 
a dependent  variable  can  be  expressed  as  an  explicit  function  of  the  other  explanatory 
variables  as  in  (1.1).  Thus,  in  the  implicit  case  the  distinction  between  response  and 
regressor  variables  is  unnecessary.  When  the  relationship  is  implicit,  then  the  parameters 
ft  cannot  be  estimated  by  most  ordinary  leeist  squares  procedures,  but  can  be  estimated 
by  orthogonal  distance  regression  procedures  such  as  ODRPACK. 

Example:  One  of  the  simplist  examples  of  an  implicit  model  is  that  of  fitting 
a circle  or  ellipse  in  a plane,  as  is  done  in  example  3.2.4  on  page  244  of  [Fuller, 
1987].  In  this  example,  shown  in  §2.C.ii,  the  data  are  observations  digitized  from 
the  x-ray  image  of  a hip  prosthesis,  where  the  variables  x j = (vj,  hi),  i = 1, . . . , n, 
are  the  vertical  and  horizontal  distances  from  the  origin,  respectively,  and  the 
implicit  model  is  that  of  the  ellipse 

fi{xi]/3)  = Psivi  - + 2p4{vi  - Pi){hi  - P2)  + - Pi)^  w 0 

for  i = 1, . . . ,n.  □ 

Finally,  sometimes  the  observations  are  multiresponse  and  thus  must  satisfy  more 
than  one  relationship  /.  This  most  commonly  arises  when  the  observations  are  complex 
variables,  although  it  can  occur  in  other  Ceises  as  well. 

Example:  The  problem  shown  in  §2.C.iii  is  an  example  of  multiresponse  data 
that  originates  because  the  underlying  data  are  complex.  The  problem  is  de- 
scribed in  Chapter  4,  and  on  pages  280-281,  of  [Bates  and  Watts,  1988].  In  this 
case,  the  response  variable  is  the  pair  of  values  representing  the  real  and  imag- 
inary parts  of  complex- valued  impedance  measurements,  Zi,  i = l,...,n,  of  a 
polymer,  and  the  explanatory  variable,  x*,  i = 1, . . . ,n,  is  the  (real- valued)  fre- 
quency. The  function  is  explicit,  representing  the  dielectric  constant  by  a general 
model  proposed  in  [Havriliak  and  Negami,  1967], 

Pi  - P2 

(1  + 


2i  ~ + 
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for  i = where  j = v^— 1.  For  ODRPACK,  this  must  be  encoded  as 

a two-term  multiresponse  problem  with  yi  € 3?^  representing  the  pair  of  values 

The  explicit,  implicit  and  multiresponse  problem  types  are  described  in  greater  detail 
in  the  following  sections.  They  are  each  easily  solved  using  ODRPACK.  The  solution 
procedures  required  for  single  and  multiresponse  data  are  essentially  the  same,  although, 
as  noted  below,  ODRPACK  does  need  to  know  that  the  data  cire  multiresponse  rather 
than  single  response  in  order  to  find  the  correct  solution.  The  solution  procedures 
required  for  explicit  and  implicit  functions,  however,  are  different,  with  the  implicit 
relationships  generally  costing  more  to  solve.  (See  §1.B.)  Thus,  it  is  important  that 
users  be  aware  of  the  differences  between  explicit  and  implicit  models,  and  between 
single  and  multiresponse  data,  and  invoke  ODRPACK  appropriately  for  their  particular 
model  and  data  type. 

l.A.i.  Explicit  Orthogonal  Distance  Regression 

We  define  the  explicit  orthogonal  distance  regression  problem  as  follows.  Let  (x^,  y^),  i = 
1, . . . , n,  be  an  observed  set  of  data,  where,  for  simplicity,  we  assume  that  x j G 3?^  and 
yi  G 31?^.  (For  a discussion  of  the  orthogonal  distance  regression  problem  for  higher 
dimensional  data,  see  §l.A.iii.)  Suppose  that  the  values  of  yi  are  a possibly  nonlinear 
function  of  Xi  and  a set  of  unknown  parameters  /3  G 3fJ^,  but  that  both  the  Xi  and  the  y, 
contain  actual  but  unknown  errors  6*  G 3?^  and  c*  G 3?^,  respectively,  where  here,  and 
in  the  remainder  of  this  document,  we  use  a superscript  to  denote  such  actual  but 
unknown  quantities.  Then  the  observed  value  of  yi  satisfies 

yi  = ^i\P*)-  ^i  i = l,...,n, 

for  some  actual  but  again  unknown  value  /?*. 

The  explicit  orthogonal  distance  regression  problem  is  to  approximate  (5*  hy  finding  the 
P for  which  the  sum  of  the  squares  of  the  n orthogonal  distances  from  the  curve  /(x;  /3) 
to  the  n data  points  is  minimized.  This  is  accomplished  by  the  minimization  problem 

+ (1.2) 

subject  to  the  constraints 

2/i  = /i(®*  + - Ci  i = l,...,n.  (1.3) 

Since  the  constraints  (1.3)  are  linear  in  e^,  we  can  eliminate  them  and  thus  the  €i  from 
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(1.2),  thereby  obtaining 


p,6  V 

We  then  generalize  (1.4)  to  the  weighted  orthogonal  distance  regression  (ODR)  problem 

n 

min  - yi\^  + (1-5) 

t=l 

by  introducing  the  weights  G 3?^  and  G 3?^,  t = 1, . . . ,n,  which  are  sets  of  non- 
negative numbers.  The  weights  and  can  thus  be  used  to  compensate  for  instances 
when  the  yi  and  Xi  have  unequal  precision,  or  when  there  are  observations  that  should 
be  eliminated  from  the  analysis.  (See  §1.F.) 

l.A.ii.  Implicit  Orthogonal  Distance  Regression 

The  univariate  implicit  orthogonal  distance  regression  problem  arises  from  the  assump- 
tion that  there  is  no  distinguished  variable  y and  thus  that  the  data  must  satisfy 

fi{xi  + 6*]P*)  = Q t=l,...,n. 

For  G 3?^ , this  yields  the  minimization  problem 

n 

i=l 


(1.6) 


subject  to  fi{xi  -f  )9)  = 0 i = 1, . . . , n, 

where  ws^  G , i = 1, . . . , n,  again  denotes  a set  of  positive  numbers  used  for  weighting. 
(See  §1.F.)  The  constraints  in  (1.6)  are  not  assumed  to  be  linear  in  6i,  and  therefore 
Ccinnot  be  eliminated  to  create  an  unconstrained  problem  as  was  done  above  to  specify 
the  explicit  problem  (1.4).  The  implicit  orthogonal  distance  regression  problem  must 
therefore  be  solved  differently  than  the  explicit  problem.  (See  §1.B.) 

l.A.iii.  Orthogonal  Distance  Regression  for  Multidimensional  Data 

For  simplicity,  the  preceding  two  sections  developed  the  explicit  and  implicit  orthogo- 
nal distance  regression  problem  for  univariate  data,  i.e.,  for  i»  G 3f?^  and  yi  G with 
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fi  : ^ 3fJ'.  In  this  section,  we  generalize  these  problems  to  allow  for  multivari- 

ate explanatory  data  x,  G J?"*,  m > 1,  and  multiresponse  dependent  variables  yi  G 5?^, 
q > 1,  which  are  modeled  by  fi  : — ♦ 5^.  For  such  higher  dimensional  data,  the 

errors  G 3?^  associated  with  x*,  and  the  errors  c*  G 3?^  associated  with  yi,  can  be 
weighted  by  wsi  € and  G 31?^^^^,  i = 1, , ..,n,  respectively,  where  the  ty^^must 

be  positive  definite  matrices  and  the  We^-must  be  positive  semidefinite  matrices.^ 

For  the  explicit  problem,  we  assume  that  the  q responses  of  the  observed  values  of  yi 
satisfies 

yi  = fiixi  + 6l;/3*)- t=l,...,n, 

for  some  actual,  but  unknown  value  /3*.  Thus,  if  we  let  X G 3?”’^*"  and  Y G denote 
the  arrays  with  tth  row  Xj  and  yi,  respectively,  and  A G and  E G 3?”^®  denote 

the  arrays  with  ith  row  Si  and  €i,  respectively,  we  cLre  then  assuming  that  the  observed 
values  X and  Y satisfy 


X = X*-A* 

Y = Y*-E*, 

and  that  the  estimated  values,  which  we  denote  by  a “hat,”  will  satisfy 

X = Z + A 

Y = Y+E. 

Generalizing  the  problem  presented  in  §l.A.i,  the  explicit  multiresponse  orthogonal  dis- 
tance regression  problem  for  multivariate  explanatory  data  is  defined  as 

.—1 


subject  to 


Yii  = fiiixi  + 6i-,P)-Eii  ' 
Yi2  = fi2{Xi  + bi\fi)  — Ei2 


i=  1, 


• ) 


n , 


Yiq  — fiq{xi  + bi]  0)  — Eiq 


(1.7) 


where  subscript  il  denotes  the  (t,  /)th  element  of  the  corresponding  array,  e.g.,  Yu  denotes 
the  lib.  response  of  the  tth  observation,  and  /ti(x,-  -|-  6i;P)  denotes  the  model  for  the  /th 

^The  matrix  A is  positive  definite ]£ nud  only  ii  = A,  and  a^Aa  > 0 for  all  nonzero  vectors  a.  It  is 
positive  semidefinite  if  and  only  if  A’’^  = A,  and  a^Aa  > 0 for  all  vectors  a.  When  A is  positive  definite 
or  positive  semidefinite,  then  there  is  an  upper  triangular  matrix  U such  that  A = if^U.  ODRPACK 
checks  for  positive  (semi)definiteness  by  attempting  to  £3rm  V. 
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response  of  observation  t.  Because  the  constraints  in  (1.7)  are  linear  in  E, 
equivalent  to  the  problem 

n 


min 


E 


([/i(®i  + - yiV'^lMxi  + Si;P)  - Vi]  + SJwsiSi) 


imnS{/3,6)  , 

p,6 


this  is 

(1.8) 


where  we  define  the  weighted  sum  of  squares  S{P,6)  by 


S{P,6)  = S.{P,6)  + S6iP,S) 


with 

S^{P,  ^)  = ([/»(xi  + Si;  P)  - yiVweilfiixi  + Si;  P)  - yi]) 

*=i 

and 

»=1 

The  implicit  multiresponse  orthogonal  distance  regression  problem  for  multivariate  ex- 
planatory data  is 

= rmnSsiPJ) 

p,6  . P,0 


fii{xi  + Si;P)  = 0 ] 

fi2ixi  + Si;P)  = 0 I 

subject  to  , > 

fiqixi-\-Si;P)  = 0 J 

As  was  the  case  in  §l.A.ii,  the  constraints  in  (1.9)  are  not  assumed  to  be  linear  in  Si, 
and  therefore  cannot  be  eliminated  to  create  an  unconstrained  problem  as  was  done 
above  to  specify  the  explicit  problem  (1.8).  The  implicit  orthogonal  distance  regression 
problem  must  therefore  be  solved  differently  than  the  explicit  problem.  (See  §1.B.) 

Note  that  when  g > 1,  the  responses  of  a multiresponse  orthogonal  distance  regression 
problem  cannot  simply  be  treated  as  q separate  observations  as  can  be  done  for  ordinary 
least  squares  when  the  q responses  are  uncorrelated.  This  is  because  ODRPACK  would 
then  treat  the  variables  associated  with  these  q observations  as  unrelated,  and  thus  not 
constrain  the  errors  Si  in  i*  to  be  the  same  for  each  of  the  q occurrences  of  the  tth 
observation.  The  user  must  therefore  indicate  to  ODRPACK  when  the  observations 
are  multiresponse,  so  that  ODRPACK  can  make  the  appropriate  adjustments  to  the 
estimation  procedure.  (See  §2.B.ii,  subroutine  argument  NQ.) 
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l.B.  Algorithm 

The  algorithm  implemented  in  ODRPACK  is  described  in  [Boggs  et  al.,  1987  and 
1989].  Briefly,  the  solution  is  found  iteratively  using  a trust  region  Levenberg-Marquardt 
method,  with  scaling  used  to  accommodate  problems  in  which  estimated  values  have 
widely  varying  magnitudes.  The  Jacobian  matrices,  i.e.,  the  matrices  of  first  partial 
derivatives  of  fui^i  + » = l,...,n,  and  I — l,...,q,  with  respect  to  each  com- 

ponent of  P and  A,  are  computed  at  every  iteration  either  by  finite  differences  or  by  a 
user  supplied  subroutine,  as  specified  by  subroutine  argument  JOB  (see  §2.B.ii,  and  also 
§4.A).  The  iterations  are  stopped  when  any  one  of  three  stopping  criteria  are  met.  Two 
of  these  indicate  the  iterations  have  converged  to  a solution.  These  are  sum  of  squares 
convergence^  which  indicates  that  the  change  in  the  weighted  sum  of  the  squared  ob- 
servation errors  is  sufficiently  small,  and  parameter  convergence,  which  indicates  the 
change  in  the  estimated  values  of  P and  A is  sufficiently  small.  The  third  stopping 
criterion  is  a limit  on  the  number  of  iterations. 

ODRPACK  finds  the  solution  of  an  implicit  orthogonal  distance  regression  problem 
using  the  classic  quadratic  penalty  function  method.  The  penalty  function  is 

= + + + . (i.io) 

1=1 

with  penalty  term 

i=\ 

and  penalty  parameter  rk-  A sequence  of  unconstrained  minimization  problems 

imnP{P,6\rk) 

p,6 

is  then  solved  for  a sequence  of  values  of  the  penalty  parameter  r*  tending  to  oo. 
These  problems  are  equivalent  to  an  explicit  orthogonal  distance  regression  problem 
with  £{  = fi{xi  -1-  6i]P)  . As  r*.  — » 00,  — » 0,  t = 1,.  ..,n,  and  the  solution  to  (1.10) 
will  approach  that  of  the  implicit  problem  defined  by  (1.9).  See,  e.g.,  [Gill  et  al.,  1981] 
for  further  discussion  of  penalty  function  methods. 

l.C.  Specifying  the  Task 

The  user  has  the  option  of  specifying  several  different  aspects  of  the  problem: 

1.  whether  the  fit  is  to  be  by  explicit  or  implicit  orthogonal  distance  regression,  or 
by  ordinary  least  squares  (see  §1.A); 
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2.  whether  the  necessary  Jacobian  matrices  should  be  approximated  by  ODRPACK 
using  forward  or  central  finite  differences,  or  whether  the  user  has  supplied  the 
code  to  compute  them,  and,  if  such  code  has  been  provided  by  the  user,  whether 
it  should  be  checked  (see  §4.A); 

3.  whether  the  covariance  matrix  Vp  and  standard  deviations  should  be  com- 
puted for  the  estimators  of  and  whether  the  Jacobian  matrices  should  be 
recalculated  at  the  solution  for  this  computation  (see  §4.B); 

4.  whether  the  errors  A have  been  initialized  by  the  user  (see  §1.E,  and  §2,B.ii, 
subroutine  arguments  JOB  and  WORK); 

5.  whether  the  fit  is  a “restart,”  i.e.,  whether  the  fit  will  use  information  preserved 
in  the  vectors  BETA,  WORK  and  IWORK  to  continue  from  a previously  found  point 
(see  Chapter  3); 

6.  whether  the  data  are  multiresponse  or  not; 

7.  whether  subsets  of  the  unknowns  0 and  A should  be  treated  as  constants  with 
their  values  held  “fixed,”  allowing  the  user  to  examine  the  results  obtained  by 
estimating  subsets  of  the  parameters  of  a general  model  without  rewriting  the 
model  subroutine,  and  allowing  the  user  to  specify  that  some  components  of  X 
are  to  be  treated  as  if  they  are  known  exactly; 

8.  whether  weighted  or  unweighted  analysis  should  be  performed  (see  §1.F);  and 

9.  whether  the  unknowns  0 and  A should  be  scaled  to  compensate  for  cases  where 
their  values  vary  widely  in  magnitude  (see  §4.D). 

The  first  5 of  these  options  are  specified  by  ODRPACK  subroutine  argument  JOB;  mul- 
tiresponse data  are  indicated  by  subroutine  argument  NQ;  parameter  “fixing”  is  specified 
by  subroutine  arguments  IFIXB  and  IFIXX;  weighting  is  controlled  by  arguments  WE  and 
WD;  and  scaling  is  controlled  by  arguments  SCLB  and  SOLD.  A detailed  discussion  of  each 
of  these  subroutine  arguments  can  be  found  in  §2.B.ii. 

l.D.  ODRPACK  Generated  Results 

Results  generated  by  ODRPACK  are  returned  to  the  user  in  four  ways: 

• the  estimated  parameter  values  0 are  returned  in  subroutine  argument  BETA; 

• the  stopping  condition  is  returned  in  subroutine  argument  INFO; 

• the  computed  results  available  at  the  time  the  job  stopped  are  returned  in  the 
vectors  specified  by  subroutine  arguments  WORK  and  IWORK;  and 

• selected  results  are  summarized  in  automatically  generated  reports. 
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Arguments  BETA,  INFO,  WORK  and  IWORK  are  discussed  in  detail  in  §2.B.ii.  The  remainder 
of  this  section  describes  ODRPACK’s  automatically  generated  reports. 

ODRPACK  can  generate  two  different  types  of  reports.  The  first  is  an  error  report,  and 
the  second  is  a computation  report.  These  are  written  to  the  logical  units  specified  by 
subroutine  arguments  LUNERR  and  LUNRPT,  respectively.  The  user  can  associate  these 
units  with  a file  using  a Fortran  OPEN  statement.  If  the  user  sets  one  or  the  other  of 
these  units  to  a negative  value,  then  the  corresponding  report  will  be  generated  on  unit 
6.  On  most  systems,  unit  6 is  “standard  output,”  which  for  an  “interactive”  run  will 
be  the  screen.  If  the  user  sets  one  or  the  other  of  these  units  to  0,  then  generation  of 
the  corresponding  report  is  suppressed. 

ODRPACK’s  computation  reports  (see  §l.D.ii  below)  can,  at  the  user’s  option,  be  writ- 
ten simiiltaneously  to  both  the  file  associated  with  the  unit  specified  by  LUNRPT  and  to 
unit  6 (assuming  LUNRPT  does  not  specify  unit  6).  This  option,  invoked  by  ODRPACK 
subroutine  argument  IPRINT,  enables  the  user  to  monitor  ODRPACK’s  progress  inter- 
actively while  still  preserving  the  results  in  a file  for  future  reference.  The  option  is 
especially  useful  when  the  user’s  model  is  time  consuming  and  the  user  is  interactively 
experimenting  with  starting  values  and  other  ODRPACK  options:  results  printed  to 
the  screen  can  be  used  to  terminate  ineffective  or  incorrect  runs,  while  results  stored  on 
the  logical  unit  specified  by  LUNRPT  can  be  used  to  preserve  the  successful  experiments. 


l.D.i.  Error  Reports 

Error  reports  identify  incorrect  user  supplied  information  passed  to  ODRPACK  that 
prevents  the  computations  from  beginning.  For  example,  an  error  report  will  be  gen- 
erated if  the  user  specifies  the  number  of  observations  to  be  an  obviously  meaningless 
number,  such  as  a negative  value.  Error  reports  are  self  explanatory,  and  are  not  dis- 
cussed further  here.  The  information  written  in  the  error  reports  is  also  encoded  and 
returned  to  the  user  in  subroutine  argument  INFO.  (See  §2.B.ii.) 


l.D.ii.  Computation  Reports 

Computation  reports  are  divided  into  three  sections: 

a.  the  initial  report, 

b.  the  iteration  report,  and 

c.  the  final  report. 

Each  of  these  can  be  either  “short”  or  “long” . 
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l.D.ii.a.  Initial  Reports 

ODRPACK’s  initial  report  identifies  the  output  for  future  reference,  and  provides  infor- 
mation that  should  enable  the  user  to  verify  that  the  problem  was  specified  correctly. 

The  short  initial  report  is  shown  in  §2.Ci.c.  It  includes 

• the  values  n,  m,  p and  g,  the  number  n of  observations  with  nonzero  weights, 
and  the  number  p of  pcirameters  ft  actually  being  estimated; 

• the  values  specified  by  ODRPACK  subroutine  arguments  JOB,  NDIGIT,  TAUFAC, 
SSTOL,  PARTOL,  and  MAXIT,  that  indicate,  respectively,  the  requested  task,  the 
number  of  “good”  digits  in  the  user’s  function  results,  the  stopping  criteria  for 
testing  for  sum  of  squares  convergence,  the  stopping  criteria  for  testing  for  pa- 
rameter convergence,  and  the  maximum  number  of  iterations  permitted;  and 

• the  weighted  sum  of  the  squared  errors  5(/5,^),  Ss{P,S)  and  Se{0,S),  evaluated 
at  the  initial  values  of  P and  A. 

The  long  initial  report,  shown  in  §2.C.ii.c  and  §2.C.iii.c,  includes  all  the  information  in 
the  short  initial  report  and,  in  addition,  includes  a brief  summary  of  the  information 
specified  for  the  function  parameters  P,  and  also  information  about  the  user  supplied 
values  for  the  X and  Y data.  This  information  includes 

• the  starting  parameter  values  Pk,k  = 1, . . . ,p,  whether  or  not  each  parameter  is 
to  be  treated  as  “fixed”,  what  “scale”  value  will  be  used,  and,  when  the  deriva- 
tives are  constructed  using  finite  differences  the  step  size  used  in  that  calculation, 
or  when  the  user  supplied  derivatives  were  checked  what  derivatives  were  identi- 
fied as  questionable; 

• the  value  of  the  first  and  last  observation  of  each  column  of  the  explanatory 
variable  X,  and  when  the  solution  is  found  by  orthogonal  distance  regression, 

o the  starting  errors  A for  each  of  these  observations, 
o whether  or  not  these  values  will  be  held  fixed  at  their  input  values, 
o the  scale  values, 

o the  diagonal  elements  of  the  weights  wsi  associated  with  eaoh  of  these  obser- 
vations, and 

o when  the  derivatives  are  constructed  using  finite  differences  the  step  size  used 
in  that  calculation,  or  when  the  user  supplied  derivatives  were  checked  what 
derivatives  were  identified  as  questionable; 

• the  value  of  the  first  and  last  observation  of  the  response  variable  Y,  and  the 
value  of  the  corresponding  diagonal  entries  of  the  error  weights  We^. 
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l.D.ii.b.  Iteration  Reports 

The  ODRPACK  iteration  reports  enable  the  user  to  monitor  the  progress  of  the  fitting 
procedure,  where  the  user  controls  at  which  iterations  these  reports  will  be  generated. 

ODRPACK  iteration  reports  are  of  the  greatest  use  in  cases  when  ODRPACK  fails  to 
find  a satisfactory  solution.  In  cases  when  ODRPACK  does  reach  a satisfactory  solution, 
the  final  report  discussed  below  summarizes  the  most  useful  information. 

The  short  iteration  report,  shown  in  §2.C.i.c,  is  a one  line  summary  of  the  results, 
listing 

• the  iteration  number; 

• the  cumulative  number  of  function  evaluations  made  by  the  end  of  the  listed 
iteration; 

• the  weighted  sum  of  the  squared  observation  errors  S{fi,6)  at  the  current  point; 

• the  actual  relative  reduction  in  S{fi,6)  at  the  most  recently  tried  step  (used  to 
check  for  sum  of  squares  convergence); 

• the  predicted  relative  reduction  in  S{P,6)  at  the  most  recently  tried  step  (used 
to  check  for  sum  of  squares  convergence); 

• the  ratio  of  the  trust  region  radius  to  the  scaled  norm  of  (3  and  A,  which  is  an 
upper  bound  on  the  relative  change  in  the  estimated  values  possible  at  the  next 
step  (used  to  check  for  parameter  convergence);  and 

• whether  the  step  was  a Gauss-Newton  step. 

The  long  iteration  summary  lists  all  of  the  information  found  in  the  short  iteration 
summary  and,  in  addition,  includes 

• the  values  of  at  the  end  of  the  current  iteration.  (At  the  last  iteration,  the  values 
listed  will  be  those  that  produced  the  actual  and  predicted  relative  reductions 
shown  only  if  the  most  recently  tried  step  did  in  fact  make  the  fit  better.  K not, 
then  the  values  of  listed  are  those  that  produced  the  best  fit.) 

The  long  iteration  report  requires  125  columns,  and  uses  [p/3]  lines  per  iteration. 


l.D.ii.c.  Final  Reports 

The  final  report  provides  information  that  allows  the  user  to  assess  the  quality  of  the 
solution  found  by  ODRPACK. 

The  short  final  report,  shown  in  §2.C.ii.c  and  §2.C.iii.c,  includes 

• the  reason  the  computations  stopped; 

• the  number  of  iterations; 


Getting  Started 


13 


• the  number  of  function  evaluations  and,  if  the  Jacobian  was  supplied  by  the  user, 
the  number  of  Jacobian  evaluations; 

• the  rank  deficiency  of  the  model  at  the  time  the  computations  stopped; 

• the  inverse  of  the  condition  number  of  the  problem  at  the  time  the  computations 
stopped  (see  §4.C); 

• the  weighted  sum  of  the  squared  errors,  S{0,6)y  Ss{Py6)  and  Se{P,6),  evaluated 
at  the  final  values  P and  A,  and  if  the  covariance  matrix  was  computed,  the 
estimated  residual  variance  of  the  fit  and  its  associated  degrees  of  freedom;  and 

• the  final  values  /3,  whether  each  value  was  “fixed”  using  subroutine  argument 
IFIXB  or  “dropped”  by  ODRPACK  because  it  induced  rank  deficiency,  and,  if  the 
covariance  matrix  Vp  and  standard  deviations  were  computed,  the  standard 
errors  and  95%  confidence  intervals  for  the  estimators  of  p.  (See  §4.B). 

The  long  final  report  is  shown  in  §2.C.i.c.  It  includes  the  same  information  as  the  short 
final  report,  and,  in  addition,  provides  the  values  for  all  Eii,i  = 1, . . . , n,  I = 1, . . . , q, 
and  Aij,i=  1, . . . ,n,  j = 1, . . . ,m. 

l.E.  Starting  Values 

The  user  must  supply  starting  values  for  the  unknowns  being  estimated,  i.e.,  for  P and 
A.  Users  familiar  with  the  nonlinear  ordinary  least  squares  problem  are  generally  aware 
of  the  importance  of  obtaining  good  starting  values  for  p.  It  is  equally  important  here. 
Good  initial  approximations  can  significantly  decreaise  the  number  of  iterations  required 
to  find  a solution;  a poor  initial  approximation  may  even  prevent  a solution  from  being 
found.  Reasonable  initial  approximations  are  often  available  from  previous  analysis  or 
experiments.  When  good  starting  values  cire  not  readily  available,  the  user  may  have 
to  do  some  preliminary  analysis  to  obtain  them.  (See,  e.g.,  [Bates  and  Watts,  1988],  or 
[Himmelblau,  1970].) 

Users  who  do  not  provide  scale  information  are  strongly  encouraged  not  to  use  zero  as 
an  initial  approximation  for  any  of  the  parameters  P since  a zero  value  can  result  in 
incorrect  scaling.  (See  §4.D).  Setting  the  initial  approximation  to  the  largest  magnitude 
that,  for  the  user’s  problem,  is  effectively  zero,  rather  than  to  zero,  will  help  to  eliminate 
scaling  problems  and  possibly  produce  faster  convergence.  For  example,  if  Pi  represents 
change  in  a cost  measured  in  millions  of  dollars,  then  the  value  10  might  be  considered 
“effectively  zero”  and  an  initial  value  of  P^  = 10  is  preferable  to  P^  = 0. 

When  using  orthogonal  distance  regression  it  is  also  important  to  have  good  starting 
values  for  the  estimated  errors  A.  The  ODRPACK  default  is  to  initialize  A to  zero, 
which  is  the  most  obvious  initial  value.  (Note  that  zero  starting  vadues  for  A do  not 


14 


Getting  Started 


cause  scaling  problems  similar  to  those  discussed  above  for  0.)  Initializing  A to  zero, 
however,  is  equivalent  to  initially  assigning  all  of  the  errors  to  the  dependent  variable 
as  is  done  for  ordinary  least  squares.  While  this  is  quite  adequate  in  many  cases,  in 
others  it  is  not.  A plot  of  the  observed  data  and  of  the  curve  described  by  the  model 
function  for  the  initial  parameters  may  indicate  whether  or  not  zero  starting  values  for 
A are  reasonable  and  may  make  it  possible  to  visually  determine  better  starting  values 
for  some  of  the  A,j.  For  example,  if  such  a plot  shows  that  the  vertical  distance  from  a 
data  point  (zj,  y*)  to  the  curve  is  far  larger  than  the  orthogonal  distance,  then  8i  should 
probably  not  be  initialized  to  zero  for  that  point.  This  may  occur  near  an  asymptote  or 
near  a local  minimum  or  maximum.  In  such  cases,  it  is  often  appropriate  to  initialize 
6i  to  the  horizontal  distance  from  the  data  point  to  the  curve.  (See  §2.B.ii,  subroutine 
arguments  JOB  and  WORK.) 

l.F.  Weights 

Weights  can  be  used  to  eliminate  observations  from  the  analysis,  to  compensate  for 
unequal  variances  or  correlations  in  the  variables  y and  x,  or  simply  to  modify  the  effect 
of  the  various  errors  e and  8 on  the  fit. 

The  weights  6 associated  with  €i,  i = 1,.  ..,n,  must  be  positive  semidefinite 
matrices,  and  the  weights  wsi  G associated  with  8i,  i=  1, . . . ,n,  must  be  positive 

definite  matrices.^  If  is  diagonal,  then  the  errors  in  the  q responses  of  yi  are  treated 
as  uncorrelated.  Similaxly,  if  ws^  is  diagonal,  then  the  errors  in  the  m elements  of  Xi  are 
treated  as  uncorrelated.  When  the  weights  are  not  diagonal  matrices,  then  the  errors 
within  an  observation  are  treated  as  correlated.  In  all  CcLses,  however,  the  errors  between 
observations  are  treated  as  uncorrelated. 

Observations  can  be  eliminated  from  the  analysis  by  setting  the  appropriate  weight 
values  w^io  zero.  This  will  produce  the  same  results  as  an  analysis  with  the  zero- 
weighted  values  removed  from  the  data  prior  to  calling  ODRPACK.  Zero  weights  can 
be  used  to  allow  for  cases  of  multiresponse  data  where  the  number  of  responses  is  not 
constant  for  all  observations.  They  also  allow  for  easy  examination  of  the  effect  of 
outliers  and  influential  data  points. 

Outliers  can  often  be  identified  by  large  values  of  c and/or  8.  Careful  checking  of  the  data 
often  leads  to  confirmation  that  the  data  are  in  error,  and  sometimes  to  a correction. 
When  a cause  for  suspicious  data  cannot  be  found,  it  may  be  advisable  to  compare  the 
analysis  with  and  without  the  questionable  data.  Caution  is  in  order  if  the  estimates  or 
conclusions  are  highly  sensitive  to  a small  amount  of  suspicious  data.  Data  that  have 
a very  high  influence  on  a fitted  curve  may  not  result  in  large  errors,  however,  even  if 


^See  footnote  1 on  page  6. 
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they  are  in  error.  In  fact,  extremely  influential  observations  may  force  the  fitted  curve 
to  be  very  close,  leading  to  very  small  residuals.  It  is  therefore  desirable  to  identify 
influential  observations  and  to  compare  the  results  obtained  with  and  without  these 
points.  Several  methods  for  detecting  influential  observations  for  ordinary  least  squares 
are  discussed  in  [Belsley  et  ai,  1980],  [Bement  and  Williams,  1969],  [Cook,  1977],  and 
[Hoaglin  and  Welsch,  1978]. 

Using  weights  to  compensate  for  unequal  error  variances  is  not  as  straightforward  as 
using  zero  weights  to  eliminate  observations  from  the  analysis.  When  the  errors  e *,  i = 
1, . . . , n,  and  tf,-,  i=  1, . . . , n,  have  covariances  cr^.  € and  a],  G.  respectively, 

that  are  known,  then  the  weights  can  be  set  using 

vki  = cia-^  and  tqj.  = , 

i = l,...,n,  where  Ci  and  C2  are  constants  selected  so  that  the  magnitudes  of  the 
weighted  errors  and  SJvJSiSi  will  be  approximately  the  same. 

In  practice,  weights  are  often  derived  from  theory  or  obtained  from  the  data  being  fit. 
Users  must  be  aware  that  incorrectly  specified  weights  can  adversely  affect  the  results. 
(See,  e.g.,  [Boggs  and  Rogers,  1990a].)  Thus,  when  the  need  for  weights  is  suspected 
and  the  error  covariances  are  not  known,  it  is  extremely  important  to  analyze  how 
the  weights  are  affecting  the  computed  results.  See  [Fuller,  1987]  for  a more  complete 
discussion  of  weights  and  their  implications  for  the  estimation  of  the  parameters  of 
measurement  error  models,  and  also  [Bates  and  Watts,  1988]  for  a procedure  that  can 
be  used  to  estimate  the  covciriance  matrix  in  the  case  of  multiresponse  nonlinear 
ordinary  least  squares. 

l.G.  Default  Values  and  Structured  Arguments 

ODRPACK  uses  default  values  and  structured  arguments  to  simplify  the  user  interface. 
The  availability  of  default  values  in  ODRPACK  means  that  the  user  does  not  have  to 
be  concerned  with  determining  values  for  many  of  the  ODRPACK  arguments  unless 
the  problem  being  solved  requires  the  use  of  nondefault  values.  Structured  arguments, 
described  below,  can  reduce  the  cimount  of  storage  space  required  for  some  arguments, 
and  the  work  required  by  the  user  to  initialize  those  arguments. 


l.G.i.  Default  Values 

Default  values  have  been  specified  for  ODRPACK  subroutine  arguments  wherever  rea- 
sonable. These  default  values  are  invoked  when  the  user  sets  the  corresponding  argu- 
ment to  any  negative  value.  Arrays  with  default  values  are  invoked  by  setting  the  first 
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element  of  the  array  to  a negative  value,  in  which  case  only  the  first  value  of  the  array 
will  ever  be  used.  This  allows  a scalar  to  be  used  to  invoke  the  default  values  of  arrays, 
thus  saving  space  and  eliminating  the  need  to  declare  such  arrays. 

Users  are  encouraged  to  invoke  the  default  values  of  arguments  wherever  possible.  These 
values  have  been  found  to  be  reasonable  for  a wide  class  of  problems.  Fine  tuning  these 
arguments  can  then  be  done  later  if  necessary. 


l.G.ii.  Structured  Arguments 

Certain  ODRPACK  arguments  specify  attributes  of  the  individual  components  of  X, 
A,  and  Y.  These  attribute  arrays  are  frequently  either  constant  for  all  components  or 
are  constant  within  each  column  and  vary  only  between  the  columns.  The  structured 
argument  facility  allows  a scalar  to  specify  an  attribute  of  an  entire  column  or  of  the 
whole  array. 

For  example,  ODRPACK  argument  WD  specifies  the  attribute  array  of  weights 
W/c^  G whose  ith  component  specifies  the  elements  of  iqjj,  i = 1,  ...,n  (cf. 

(1.8)  ajid  (1.9)).  Suppose  X contains  temperature  measurements,  where  each  row  indi- 
cates a different  hour  at  which  the  readings  were  taken,  and  each  column  a different  day. 
Then  the  user  might  want  to  weight  each  component  of  A equally,  and  thus  would 
be  constant  throughout.  If  one  column  of  X contained  hourly  temperature  readings  and 
the  other  hourly  humidity  readings,  however,  then  the  user  would  probably  not  want  to 
weight  the  errors  in  the  temperature  measurements  the  Scime  as  those  in  the  humidity 
measurements,  but  might  want  to  specify  the  same  weight  for  each  observation.  In 
this  case,  the  components  of  would  be  constant  for  each  row  t,  i = l,...,n.  Of 
course,  in  other  cases,  the  user  might  want  to  weight  each  component  of  A differently, 
and  thus  each  component  of  would  be  different. 

ODRPACK  structured  arguments  exploit  this  structure.  If  each  of  the  elements  of  an 
attribute  array  is  the  same,  then  a scalar  can  be  used  to  specify  all  elements  of  the  array. 
If  the  values  of  such  an  array  only  vary  between  the  columns,  then  each  column  of  the 
array  can  be  specified  by  a single  value  using  a row  vector.  Thus,  it  is  only  necessary 
for  the  user  to  supply  all  elements  of  an  attribute  array  when  the  elements  of  one  or 
more  of  the  columns  must  be  individually  specified.  The  use  of  structured  arguments 
in  described  in  detail  in  §2.B.ii.  (See,  e.g.,  subroutine  argument  WD.) 
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2. A.  Subroutine  Declaration  and  Call  Statements 

The  declaration  and  call  statements  for  ODRPACK ’s  user  callable  routines  are  given 
below.  DODR  and  DODRC  invoke  the  double  precision  version  of  the  code  and  SODR  and 
SODRC  invoke  the  single  precision  version.  DODR  and  SODR  preset  many  arguments  to 
their  default  values  and  therefore  have  shorter  call  statements  than  DODRC  and  SODRC. 
In  contrast,  DODRC  and  SODRC  have  expanded  call  statements  that  give  the  user  greater 
control  in  solving  the  orthogonal  distance  regression  problem. 

The  information  in  this  section,  §2.A,  is  provided  primarily  for  reference.  Users  are 
directed  to  §2.B  for  definitions  of  the  subroutine  arguments,  and  to  §2.C  for  sample 
programs.  The  examples,  which  use  Fortran  PARAMETER  statements  to  dimension  ODR- 
PACK arrays,  provide  a recommended  format  for  creating  an  ODRPACK  driver  that 
will  allow  future  changes  to  be  made  easily. 

Note  that  although  ODRPACK  is  distributed  in  both  single  precision  and  double  pre- 
cision versions,  both  versions  may  not  be  available  to  all  users.  In  addition,  even  when 
both  versions  are  available,  the  single  precision  version  may  not  be  appropriate  to  use. 
This  is  because  ODRPACK  is  sensitive  to  the  machine’s  precision,  and  requires  approx- 
imately 14  decimal  places.  Somewhat  fewer  places  should  still  work,  but  six  or  seven 
decimal  places  are  definitely  too  few  for  general  use,  since  only  the  simplest  problems 
could  be  solved  correctly  at  such  reduced  precisions.  When  both  versions  are  available, 
the  user  must  choose  which  version  of  ODRPACK  to  use  based  upon  which  version  sup- 
plies adequate  precision  on  the  target  machine.  To  our  knowledge,  at  present  only  Cray 
and  CDC  machines  offer  sufficient  precision  to  permit  general  use  of  the  single  precision 
version  of  ODRPACK.  For  other  machines,  we  recommend  the  double  precision  version. 

If  both  versions  of  ODRPACK  have  sufficient  precision  on  the  user’s  machine,  then 
either  may  be  used.  When  both  the  single  and  double  precision  versions  are  available, 
however,  there  are  generally  trade-offs  between  them.  The  double  precision  version 
will  offer  greater  accuracy  in  results,  while  the  single  precision  version  will  require  less 
storage  and  possibly  less  machine  time. 
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DODR:  Compute  the  explicit  or  implicit  weighted  orthogonal  distance  regression,  or  linear  or 
nonlinear  ordinary  le£ist  squares  solution  in  double  precision.  Derivatives  are  either 
supplied  by  the  user  or  numerically  approximated  by  ODRPACK.  Control  variables  are 
preset  to  their  default  values,  and  a three  part  report  of  the  results  can  be  optionally 
generated. 


PROGRAM  HAIR 

I 

EXTERRAL 
+ FCR 
IRTEGER 
+ R.M.RP.RQ, 

+ LDY.LDX,  LDWE,LD2WE,LDWD,LD2WD, 

+ JOB,  IPRIRT.LURERR.LURRPT, 

+ LVORK.LIVORK,  IRFO 
IRTEGER 

IWORK(LIVORK) 

DOUBLE  PRECISIOR 

+ BETA(RP),  Y(LDY,Rq),X(LDX,M), 

+ WE(LDWE,LD2WE.RQ),tfD(LDWD,LD2WD.M). 
-I-  WORK(LWORK) 

I 

CALL  DODR 
+ (FCR, 

+ R,M,RP,RQ, 

+ BETA, 

+ Y,LDY,X,LDX, 

+ VE,LDVE,LD2WE,WD,LDVD,LD2UD, 

+ JOB, 

+ IPRIRT,LURERR,LURRPT, 

+ VORK,LVORK,IWORK,LIWOBK, 

+ IRFO) 
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DODRC:  Compute  the  explicit  or  implicit  weighted  orthogonal  distance  regression,  or  linear  or 
nonlinear  ordinary  least  squares  solution  in  double  precision.  Derivatives  are  either 
supplied  by  the  user  or  numerically  approximated  by  ODRPACK.  Control  values  are 
supplied  by  the  user,  and  a three  part  report  of  the  results  can  be  optionally  generated. 

PROGRAM  MAIH 

I 

EXTERHAL 
+ FCH 
IVTEGER 
+ I.H.IP.HQ. 

+ LDY.LDX,  LDWE,LD2WE,IDWD,LD2WD,  LDIFX, 

+ JOB.HDIGIT,  MAXIT,  IPRIHT.LUHERR.LUHRPT, 

+ LDSTPD,  LDSCLD,  LWORK , LItf ORK , IHFO 
IVTEGER 

+ IFIXB(FP),IFIXX(LDIFX,M),  IWORK(LIWORK) 

DOUBLE  PRECISIOF 
+ TAUFAC,  SSTOL.PARTOL 
DOUBLE  PRECISIOF 

+ BETA(FP),  Y(LDY,FQ),X(LDX,M). 

+ WE(LDWE,LD2WE,FQ),WD(LDWD,LD2WD,M), 

+ STPB(FP),STPD( LDSTPD, H),  SCLB(FP) ,SCLD(LDSCLD,M) . 

■«-  VORK(LVORX) 

I 

CALL  DODRC 
+ (FCF, 

+ F,M,FP,FQ, 

+ BETA, 

+ Y,LDY,X,LDX. 

VE,LDUE,LD2UE,VD,LDVD,LD2WD, 

+ IFIXB.IFIXX, LDIFX, 

+ JOB, FDIGIT, TAUFAC, 

+ SSTOL , P ARTOL , MAXIT , 

+ IPRIFT,LUFERR,LUFRPT, 

+ STPB,STPD, LDSTPD, 

+ SCLB,SCLD, LDSCLD, 

-I-  WORK,  LVORK.IWORK,  LI  WORK, 

+ IFFO) 

I 

EFD 
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SODR:  Compute  the  explicit  or  implicit  weighted  orthogonal  distance  regression,  or  linear  or 
nonlinear  ordinary  least  squares  solution  in  single  precision.  (SODR  is  appropriate  for 
general  use  only  on  mcichines  with  approximately  14  decimal  places  of  precision  for  single 
precision.)  Derivatives  are  either  supplied  by  the  user  or  numerically  approximated  by 
ODRPACK.  Control  variables  are  preset  to  their  default  values,  and  a three  part  report 
of  the  results  Ccin  be  optionally  generated. 

PROGRAM  MAIH 

I 

EXTERHAL 
+ FCI 
IHTEGER 
+ I, M, HP, HQ, 

+ LDY,LDX,  LDWE,LD2HE,LDWD,LD2WD, 

+ JOB,  IPRIHT,LUHERR,LUHRPT, 

+ LW0RK,LIH0RK,  IHFO 
IHTEGER 

-I-  IHORK(LIHORK) 

REAL 

+ BETA(HP),  Y(LDY,HQ),X(LDX,M), 

+ WE(LDWE,LD2HE,HQ),WD(LDHD,LD2HD,M), 

-I-  VORK(LVORK) 

I 

CALL  SODR 
+ (FCH, 

+ H,M,HP,HQ, 

+ BETA, 

+ Y,LDY,X,LDX, 

+ WE,LDVE,LD2WE,VD,LDHD,LD2WD, 

+ JOB, 

+ IPRIHT,LUHERR,LUHRPT, 

* V0RK,LV0RK,IV0RK,LIW0BK, 

+ IHFO) 
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SODRC:  Compute  the  explicit  or  implicit  weighted  orthogonal  distance  regression,  or  linear  or 
nonlinear  ordinary  le£ist  squares  solution  in  single  precision.  (SODRC  is  appropriate  for 
general  use  only  on  machines  with  approximately  14  decimal  places  of  precision  for  single 
precision.)  Derivatives  are  either  supplied  by  the  user  or  numerically  approximated  by 
ODRPACK.  Control  values  are  supplied  by  the  user,  and  a three  part  report  of  the 
results  can  be  optionally  generated. 

PROGRAM  HAIR 

I 

EXTERNAL 
+ FCF 
INTEGER 
+ N.M.NP.NQ, 

+ LDY.LDX,  LDWE.LD2WE.LDWD,LD2WD,  LDIFX, 

+ JOB.NDIGIT,  HAXIT,  IPRINT.LUNERR.LUNRPT, 

+ LDSTPD,  LDSCLD,  LWORK , LIWORK , INFO 
INTEGER 

+ IFIXB(NP),IFIXX(LDIFX,M),  IWORK(LIWORK) 

REAL 

+ TAUFAC,  SSTOL.PARTOL 
REAL 

+ BETA(NP),  Y(LDY,Nq),X(LDX,M), 

+ VE(LDVE.LD2WE.Nq).WD(LDVD.LD2WD.M), 

+ STPB(NP),STPD (LDSTPD, M),  SCLB(NP) ,SCLD(LDSCLD,M) , 

+ WORK (LWORK) 

I 

CALL  SODRC 
+ (FCN, 

+ N.M.NP.Nq, 

+ BETA, 

+ Y,LDY,X,LDX, 

■I-  VE,LDVE,LD2WE,VD,LDWD,LD2WD, 

+ IFIXB,IFIXX, LDIFX, 

+ JOB, NDIGIT, TAUFAC, 

+ SSTOL,PARTOL,MAXIT, 

+ IPRINT,LUNERR,LUNRPT, 

+ STPB,STPD, LDSTPD, 

+ SCLB,SCLD, LDSCLD, 

+ WORK, LWORK, IWORK, LI WORK, 

+ INFO) 


END 
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2.B.  Subroutine  Argument  Descriptions 
2.B.i.  Synopsis 

The  arguments  of  the  ODRPACK  user  callable  subroutines  are  logically  grouped  as 
shown  below.  Arguments  shown  in  parenthesis  (...)  are  not  included  in  the  DODR  and 
SODR  call  statements;  DODR  and  SODR  automatically  preset  these  variables  to  the  default 
values  given  in  §2.B.ii.  All  other  arguments  are  common  to  all  ODRPACK  user  callable 


subroutines. 

Argument 

Number  Arguments 

Group  Description 

1 FCH, 

Name  of  user  supplied  subroutine 
for  function  and  derivative  computation 

2 to  5 

Problem  size  specification 

6 BETA, 

Function  parameters 

7 to  10  Y,LDY,X,LDX, 

Dependent  and  explanatory  variables 

11  to  16  WE,LDWE,LD2WE.VD.LDWD,LD2WD, 

Weights 

17  to  19  (IFIXB.IFIXX.LDIFX,) 

Parameter  cind  variable  fixing 

20  to  22  JOB,  (IIDIGIT,TAUFAC, ) 

Computation  and  initialization  control 

23  to  25  ( SSTOL , P ARTOL , MAXIT , ) 

Stopping  criteria 

26  to  28  IPRIFT , LUHERR , LUHRPT , 

Print  control 

29  to  31  (STPB,STPD,LDSTPD,) 

Derivative  step  sizes 

32  to  34  (SCLB,SCLD,LDSCLD,) 

Scaling 

35  to  38  WORK , LVORK , IWORK , LIUORK , 

Work  vectors  and  returned  results 

39  IFFO 

Stopping  condition 
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ODRPACK  Subroutine  Argument  Definitions 

The  arguments  of  ODRPACK’s  user  callable  subroutines  are  described  below  in  order 
of  their  occurrence  in  the  call  statements.  Appropriate  declaration  statements  for  each 
argument  are  shown  in  brackets  [. . .]  following  the  argument  name;  the  character  string 
<real>  denotes  REAL  when  using  single  precision  subroutines  SODR  and  SODRC  (which 
should  be  used  only  on  machines  with  approximately  14  decimal  digits  of  precision  in 
single  precision),  and  denotes  DOUBLE  PRECISION  when  using  double  precision  subrou- 
tines DODR  and  DODRC.  Each  argument  is  numbered  as  shown  in  §2.Bi,  and  will  be  cross 
referenced  by  both  number  and  name,  e.g.,  1-FCN,  enabling  the  user  to  easily  find  the 
definition  of  a specific  argument.  In  axidition,  a flag  indicating  whether  the  argument 
is  passed  to  or  from  ODRPACK  is  placed  in  the  left  margin  by  the  argument  number. 
The  flags  are: 

indicating  the  eirgument  specifies  information  that  must  be  input  to  ODRPACK, 
and  that  the  input  information  is  preserved  upon  return  from  the  ODRPACK 
subroutine; 

•<=  indicating  the  argument  specifies  information  that  is  output  by  the  ODRPACK 
subroutine;  and 

indicating  the  argument  specifies  information  that  must  be  input  to  ODRPACK, 
but  that  the  input  values  will  be  overwritten  by  ODRPACK  upon  return  from 
the  subroutine. 


NOTE 

Substitute  DOUBLE  PRECISION  for  <real>  when  using  DODR  and  DODRC. 
Substitute  REAL  for  <real>  when  using  SODR  and  SODRC. 
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=^1  - FCN  [EXTERNAL  FCN] 

The  name  of  the  user  supplied  subroutine  that,  given  the  current  values  of  the 
explanatory  variable,  Xi+^i,  and  the  current  values  of  the  function  parameters, 
/S,  computes  the  predicted  values 


F(I,L)  = /il(xi  + ^i;/3), 

I = 1, . . . , 71,  and  L = 1, . . . , g,  and  optionally  computes  the  matrices  of  first 
partial  derivatives,  i.e.,  the  Jacobian  matrices 


FJACB(I,K,L) 


9fn,{^i  ~l~  ^iiP) 


for  I = 1, ...  ,71,  K = 1, . . .,p,  and  L = 1, . . .,g,  and 


FJACD(I,J,L) 


dfri,{xi  + Si]P) 


for  I = 1, . . . , 71,  J = 1, . . . , 771,  and  L = 1, . . . , g.  The  code  for  evaluating  F 
must  always  be  provided.  The  code  for  evaluating  FJACB  and  FJACD  is  required 
only  when  the  second  digit  of  argument  20- JOB  is  greater  than  or  equal  to 
two.  (When  the  second  digit  of  JOB  is  zero  or  one,  the  necessary  Jacobian 
matrices  will  be  approximated  by  ODRPACK  using  finite  differences.  See 
argument  20-JOB,  and  also  §4.A.)  The  code  for  FJACD  does  not  need  to  be 
supplied  when  the  fit  is  by  OLS. 

At  a given  call  to  subroutine  FCN,  ODRPACK  will  never  request  that  both  the 
function  values  (F)  and  the  derivative  values  (FJACB  and  FJACD)  be  computed. 
While  it  is  generally  most  cost  effective  if  the  user  only  performs  the  required 
computations,  it  is  not  an  error  for  both  function  values  and  derivatives  to  be 
computed  each  time  FCN  is  invoked.  Note,  however,  that  array  FJACD  must 
never  he  altered  when  the  fit  is  by  ordinary  least  squares,  since  no  space  is 
assigned  to  that  array  in  that  case. 

The  user  must  supply  a value  for  every  element  of  the  selected  arrays.  If  some 
responses  of  some  observations  are  actually  missing,  then  the  user  can  set  the 
corresponding  weights  in  argument  11-WE  to  zero  in  order  to  remove  the  effect 
of  the  missing  observation  from  the  analysis.  (See  §1.F.) 

ODRPACK ’s  parameter  fixing  arguments  IFIXB,  IFIXX  and  LDIFX  are  passed 
to  subroutine  FCN  for  the  user’s  convenience.  When  a parameter  is  fixed,  then 
the  derivative  with  respect  to  that  parameter  can  be  set  to  zero.  ODRPACK 
will  automatically  zero  these  derivative  values,  however,  and  thus  it  is  not 
necessary  for  the  user  to  be  concerned  with  this  unless  the  derivative  compu- 
tations are  especially  expensive.  (See  Chapter  3.) 
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The  argument  list  cind  dimension  statements  for  subroutine  FCN  must  be  ex- 
actly as  shown  below. 


SUBROUTIFE  FCH(H,M,HP,FQ, 


+ 


+ 


+ 


+ 


+ 


LDI.LDM.LDIP. 
BETA.XPLUSD, 
IFIXB.IFIXX.LDIFX, 
IDEVAL , F , F JACB , FJACD , 
ISTOP) 


C 

C 


C 


C 


C 


C 


ISPUT  ARGUMENTS 

(WHICH  MUST  HOT  BE  CHANGED  BY  THIS  ROUTINE) 

INTEGER  IDEVAL , LDFIX , LDM , LDN , LDNP , M , N , NP , HQ 
INTEGER  IFIXB(NP),IFIXX(LDIFX.M), 

<real>  BETA(NP) ,XPLUSD(LDN,M) 

OUTPUT  ARGUMENTS 
INTEGER  ISTOP 

<real>  F(LDN,Nq) , FJACB (LDN, LDNP, NQ) , FJACD (LDN, LDM, NQ) 

< set  ISTOP  > 

IF  (ISTOP. HE. 0)  RETURN 

COMPUTE  FUNCTION 

IF  (MOD (IDEVAL, 10 ).GE.l)  THEN 

< compute  F(I,L),  1=1,..., N,  A L=1,...,NQ  > 

END  IF 

COMPUTE  DERIVATIVES  WITH  RESPECT  TO  BETA 
IF  (M0D(IDEVAL/10,10).GE.l)  THEN 

< compute  FJACB(I,K,L),  1=1,..., H,  K=1,...,NP,  A L=1,...,NQ  > 
END  IF 

COMPUTE  DERIVATIVES  WITH  RESPECT  TO  DELTA 
IF  (M0D(IDEVAL/100,10) .GE.l)  THEN 

< compute  FJACD(I,J,L),  1=1,..., N,  J=1,...,M,  A L=l,...,Nq  > 


C 

C 


C 


C 


C 


C 


END  IF 


RETURN 

END 


where 
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N is  the  number  of  observations,  n. 

M is  the  number  of  elements,  m,  in  each  explanatory  variable,  x j 6 3?”*, 
i.e.,  the  number  of  columns  of  data  in  X + A. 

NP  is  the  number  of  function  pzirameters,  p. 

NQ  is  the  number  of  responses,  q,  per  observation. 

LDN  is  an  array  leading  dimension  declarator  that  equals  or  exceeds  n. 

LDM  is  an  array  leading  dimension  declarator  that  equals  or  exceeds  m. 

LDNP  is  an  array  leading  dimension  declarator  that  equals  or  exceeds  p. 

BETA  is  the  singly  subscripted  array  that  contains  the  current  values  of  /3. 

XPLUSD  is  the  doubly  subscripted  array  that  contains  the  current  value  of 
the  explanatory  variables,  i.e.,  X + A. 

IFIXB  is  the  singly  subscripted  array  designating  whether  a given  parameter 
Pk  is  to  be  treated  as  fixed.  (See  argument  17-IFIXB  below.) 

IFIXX  is  the  doubly  subscripted  array  designating  whether,  when  the  solu- 
tion is  found  by  orthogonal  distance  regression,  a given  explanatory 
variable  Xu  is  to  be  treated  as  without  error.  (See  argument  18- 
IFIXX  below.) 

LDIFX  is  the  leading  dimension  of  array  IFIXX.  (See  argument  19-LDIFX 
below.) 

IDEVAL  is  a three  digit  INTEGER  variable  with  decimal  expansion  J3J2J1 
passed  to  subroutine  FCN  by  ODRPACK  to  designate  what  values 
are  to  be  computed. 


Computation 

Digit 

F 

Ji  = 0 array  F need  not  be  computed 

= 1 function  values  F must  be  computed 

(for  constructing  finite  difference  derivatives) 
= 2 function  values  F must  be  computed 
(for  evaluating  S{0,6)  at  new  point) 

= 3 function  values  F must  be  computed 
(for  miscellaneous  calculations) 

FJACB 

J2  = 0 array  FJACB  need  not  be  computed 

= 1 derivative  values  FJACB  must  be  computed 

FJACD 

J3  = 0 array  FJACD  must  not  be  altered 

= 1 derivative  values  FJACD  must  be  computed 
(required  for  ODR  fits  only) 

F is  the  doubly  subscripted  array  in  which  the  nx  q array  of  predicted 
values  for  each  response  of  each  observation  must  be  stored. 
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F JACB  is  the  triply  subscripted  array  in  which  the  nxpx  q array  of  partial 
derivatives  with  respect  to  P for  each  response  of  each  observation 
must  be  stored. 

F JACD  is  the  triply  subscripted  array  in  which  the  nxmxq  array  of  partial 
derivatives  with  respect  to  A for  each  response  of  each  observation 
must  be  stored. 

ISTOP  is  a variable  that  enables  the  user  to  instruct  ODRPACK  to  reject 
the  current  values  in  BETA  cind  XPLUSD  as  unacceptable. 


ISTOP 

Meaning 

< 0 
= 0 

> 0 

regression  procedure  should  be  stopped 
current  ^ and  A are  acceptable  for  use  by  FCH: 

requested  values  were  properly  computed  within  FCI  and 
regression  procedure  should  continue 
current  and  A are  unacceptable  for  use  by  FCS: 
if  call  to  FCI  was  for  constructing  derivatives 
regression  procedure  will  be  stopped 
if  call  to  FCI  was  for  evaluating  5(0,6) 

a new  point  will  be  selected  closer  to  the  most 
recently  tried  acceptable  values  of  0 and  A 
if  call  to  FCI  was  for  miscellaneous  calculations 
regression  procedure  will  be  stopped 

Note  that  even  when  ISTOP  is  used  to  stop  the  regression  procedure, 
the  final  summary  of  the  computation  report  will  still  be  generated 
if  it  ha.s  been  requested.  (See  argument  26-IPRINT.) 


=►2  - N [INTEGER  N] 

The  number  of  observations,  n. 

=>3  - M [INTEGER  M] 

The  ntunber  of  elements,  m,  in  each  explanatory  variable  € 3?”^,  i.e.,  the 
number  of  columns  of  data  in  X. 

=►4  - NP  [INTEGER  NP] 

The  number  of  function  parameters,  p. 

==J^5  - NQ  [INTEGER  NQ] 

The  number  of  responses,  q,  per  observation. 

^6  - BETA  [<real>  BETA(NP)] 

The  singly  subscripted  array  that  contains  the  (current)  values  of  /3. 
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On  input:  BETA  must  contain  initial  approximations  for  the  function  param- 
eters p.  Initial  approximations  should  be  chosen  with  care  since 
poor  initial  approximations  can  significantly  increase  the  number 
of  iterations  required  to  hud  a solution  and  possibly  prevent  the 
solution  from  being  found  at  all.  (See  §1.E.) 

On  return:  BETA  contains  the  “best”  estimate  of  the  parameters,  at  the 
time  the  computations  stopped. 

=^7  -Y  [<real>  Y(LDY,NQ)] 

The  double  subscripted  array  that  contains  the  values  of  the  response  variable 
Yu.,1  = 1,...  ,n,  L = l,...,g.  When  the  model  is  explicit,  the  user  must 
supply  a value  for  each  of  the  n x q elements  of  Y;  if  some  responses  of  some 
observations  are  zu:tually  missing,  then  the  user  can  set  the  corresponding 
weight  in  argument  11-WE  to  zero  in  order  to  remove  the  effect  of  the  missing 
observation  from  the  analysis.  When  the  model  is  implicit,  Y is  not  referenced. 
(See  argument  20-JOB,  and  §1.F.) 

=>8  - LDY  [INTEGER  LDY] 

The  leading  dimension  of  array  Y.  AVhen  the  model  is  explicit,  LDY  must  equal 
or  exceed  n.  When  the  model  is  implicit,  LDY  must  equal  or  exceed  1. 

=>9  - X [<real>  X(LDX,M)] 

The  doubly  subscripted  array  that  contains  the  observed  values  of  the  ex- 
planatory variable  X. 

=>10  - LDX  [INTEGER  LDX] 

The  leading  dimension  of  array  X.  LDX  must  equal  or  exceed  n. 

=^>11  - WE  [<real>  WE(LDWE,LD2WE,Nq)] 

The  triply  subscripted  array  that,  when  the  model  is  explicit  specifies  how 
each  €i  is  to  be  weighted  in  the  weighted  orthogonal  distance,  and  when  the 
model  is  implicit  specifies  the  starting  penalty  parameter  value,  tq  (see  §1.A, 
§1.B,  §1.F,  and  argument  20-JOB). 

For  explicit  models,  WE  is  a structured  argument:  only  the  specific  elements 
of  WE  identified  in  the  table  below  are  referenced  by  ODRPACK.  (See  §1.G.) 
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WE(1,1,1)  LDWE  LD2WE 

For  I = 1,. . .,n, 

t«ei 

<0  — — 

= — WE(1,1,1)J,,  Ig  a,  q X q identity  matrix 

>0  =1  =1 

= diag{WE(l.l,L2),  L2  = 1,...,9} 

> n =1 

= diag{WE(I,l.L2),  L2  = 1,...,5} 

= 1 >9 

= WE(1,L1,L2),  LI  = l,...,g,  & L2  = 1,...,9 

> n > q 

= WE(I,L1,L2),  LI  = 1,...,9,  & L2  = 

For  implicit  models,  only  the  first  element  of  WE  is  ever  referenced,  and  tq  is 
set  as  follows.  


WECl.l.l) 

ro 

o o 

VI  A 

= 10 

= WECl.l.l) 

=J^12  - LDWE  [INTEGER  LDWE] 

The  leading  dimension  of  array  WE.  LDWE  must  either  exactly  equal  one,  or 
must  equal  or  exceed  n.  See  argument  11-WE  for  further  details. 

=>13  - LD2WE  [INTEGER  LD2WE] 

The  second  dimension  of  array  WE.  LD2WE  must  either  exactly  equal  one,  or 
must  equal  or  exceed  q.  See  argument  11-WE  for  further  details. 

=>14  - WD  [<real>  WD(LDWD,LD2WD,M)] 

The  triply  subscripted  array  that  specifies  how  eaoh  6i  is  to  be  weighted  in 
the  weighted  orthogonal  distance.  (See  §1.A  and  §1.F.)  WD  is  a structured 
argument;  only  the  specific  elements  of  WD  identified  in  the  table  below  are 
ever  referenced  by  ODRPACK.  (See  §1.G.) 


WD(l.l.l)  LDWD  LD2WD 

For  I = 1, . . ., n, 
ws, 

<0  — — 

= — WD(l.l.l)/,ni  Jm  Au  X m identity  matrix 

= 0 — — 

= Imi  Im  an  m X m identity  matrix 

>0  =1  =1 

= diag{WD(l.l.J2),  J2  = l,...,m} 

> n =1 

= diag{WD(I.l.J2),  J2  = l,...,Tn} 

= 1 > m 

= WD(1.J1.J2),  J1  = l,...,m,&J2  = l,...,m 

> n ">  m 

= WD(I.J1.J2),  J1  = l,...,m,&J2  = l,...,Tn 

=»15  - LDWD  [INTEGER  LDWD] 

The  leading  dimension  of  array  WD.  LDWD  must  either  exeictly  equal  one,  or 
must  equal  or  exceed  n.  See  argument  14-WD  for  further  details. 
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=j>16  - LD2WD  [INTEGER  LD2WD] 

The  second  dimension  of  array  WD.  LD2WD  must  either  exactly  equal  one,  or 
must  equal  or  exceed  m.  See  argument  14-WD  for  further  details. 

=>17  - IFIXB  [INTEGER  IFIXB(NP)] 

The  singly  subscripted  array  that  specifies  the  indicator  variable  des- 
ignating whether  /3k  is  to  be  treated  as  “fixed,”  i.e.,  is  to  be  treated  as  a 
constant,  or  is  to  be  “unfixed”  and  thus  is  to  be  estimated. 

• If  — 0 then  /3k  is  fixed  and  BETA(K)  is  not  changed. 

• If  ^ 0 then  /3k  is  unfixed  cind  BETA(K)  is  overwritten  by 

The  default  values  are  = 1,  K = 1, . . .,p. 


IFIXB (1) 

ForK= 

< 0 
> 0 

= 1 

= IFIXB (K) 

=>18  - IFIXX  [INTEGER  IFIXX(LDIFX,M)] 

The  doubly  subscripted  array  that  specifies  the  indicator  variable  desig- 
nating whether,  when  the  solution  is  found  by  orthogonal  distance  regression, 
the  (I,  J)th  element  of  the  explanatory  variable  Xu  is  to  be  treated  as  with- 
out error  and  thus  Au  is  to  be  “fixed”  at  zero,  or  whether  that  observation  is 
“unfixed”  and  therefore  the  error  Au  is  to  be  estimated.  (When  the  solution 
is  foimd  by  ordinary  least  squares,  the  Xu  we  always  treated  as  fixed,  and 
thus  Au  = 0 for  all  I = 1, . . . , n,  and  J = 1, . . . , m.) 

• K = 0 then  Xu  is  fixed  and  Au  is  set  to  zero. 

• K ^ 0 then  Xu  is  unfixed  and  Au  is  estimated. 

The  default  values  are  = 1,  I = 1, . . . , n,  and  J = 1, . . . , m. 

IFIXX  is  a structured  argument:  only  the  specific  elements  of  IFIXX  identified 
in  the  table  below  are  referenced  by  ODRPACK.  (See  §1.G.) 


IFIXX (1,1)  LDIFX 

For  I = l,...,n, 

& J = l,...,m, 

<0  — 

>0  =1 

> n 

= 1 

= IFIXXd.J) 

= IFIXXd.J) 

=^19  - LDIFX  [INTEGER  LDIFX] 

The  leading  dimension  of  array  IFIXX.  LDIFX  must  either  exactly  equal  one. 
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or  must  equal  or  exceed  n.  See  argument  18-IFIXX  for  further  details. 

=>20  - JOB  [INTEGER  JOB] 

The  variable  that  specifies  problem  initialization  and  computational  methods. 
(See  §1.C.)  The  default  options,  selected  when  JOB  < 0,  are  that 

• the  solution  will  be  found  by  explicit  orthogonal  distance  regression; 

• the  derivatives  will  be  computed  by  forward  finite  differences; 

• the  covariance  matrix  will  be  computed  using  Jacobian  matrices  recalcu- 
lated at  the  solution; 

• A will  be  initialized  to  zero;  and 

• the  fit  will  not  be  a restart. 

When  JOB  > 0,  it  is  assumed  to  be  a 5 digit  INTEGER  with  decimal  expansion 
J5J4Z3J2J1,  where  each  digit  controls  a different  option. 


Option 

Digit  Selection 

Computational  method 
(see  §1.A) 

Ji  = 0 explicit  orthogonal  distance  regression 
= 1 implicit  orthogonal  distance  regression 
> 2 ordinary  least  squares 

Derivative  calculation 
(see  §4.A) 

I2  = 0 forward  Unite  differences 
= 1 central  finite  differences 
= 2 user  supplied  derivative  code, 
checked  by  ODRPACK 
> 3 user  supplied  derivative  code, 
not  checked  by  ODRPACK 

Covariance  matrix  Vp, 
iz.  standard  deviation  ap 
(see  §4.B) 

X3  = 0 Vp  and  ap  calculated  using 

derivatives  recomputed  at  solution 
= 1 Vp  and  <Tp  calculated  using 

derivatives  from  last  iteration 
>2  Vp  and  ap  not  calculated 

A Initialization 
(see  §1.E) 

I4  = 0 A initialized  to  zero  by  ODRPACK 
>1  A initialized  by  user 

(see  argument  36-WORK) 

Restart  facility 
(see  Chapter  3) 

J5  = 0 fit  is  not  a restart 
>1  fit  is  a restart 

=>21  - NDIGIT  [INTEGER  NDIGIT] 

The  variable  that  specifies  the  number  of  reliable  decimal  digits  tj}  in  the 
values  computed  using  subroutine  FCN.  The  value  ^ is  needed  to  calculate  the 
default  values  for  the  relative  step  sizes  used  in  calculating  finite  difference 
derivatives.  (See  arguments  29-STPB  and  30-STPD,  and  §4.A.)  It  is  also  used 
to  determine  when  the  Jacobian  with  respect  to  one  or  more  of  the  parameters 


32 


Using  ODRPACK 


P appears  to  be  rank  deficient.  It  can  not  exceed  the  number  of  decimal 
digits  $ carried  by  the  user’s  computer  for  a <real>  value.  The  default  value 
is  experimentally  determined  by  ODRPACK  for  the  user’s  particular  model. 
This  determination  of  rj;  requires  4 evaluations  of  the  model  function  from 
user  supplied  subroutine  FCN. 


IDIGIT 

< 1 
> 2 

= default  value 
= min{HDIGIT, 

=>22  - TAUFAC  [<real>  TAUFAC] 

The  variable  that  specifies  the  factor  r used  to  initialize  the  trust  region 
radius.  The  trust  region  is  the  region  in  which  the  local  approximation  to 
S{P,  6)  is  considered  reliable.  The  diameter  of  this  region  is  adaptively  chosen 
at  each  iteration  based  on  information  from  the  previous  iteration.  At  the  first 
iteration,  the  initial  diameter  is  set  to  r times  the  length  of  the  full  Gauss- 
Newton  step  calculated  at  the  initial  values  of  ft  and  A.  The  default  value  is 
r = 1.  When  r < 1 the  size  of  the  initial  step  attempted  at  the  first  iteration 
is  smaller  than  the  full  Gauss-Newton  step.  This  may  be  appropriate  if,  at 
the  first  iteration,  the  computed  results  for  the  full  Gauss-Newton  step  cause 
an  overflow,  or  cause  and/or  A leave  the  region  of  interest. 


TAUFAC 

T 

< 0 
>0 

= 1 

= mm{TAUFAC,  1} 

=>23  - SSTOL  [<real>  SSTOL] 

The  variable  that  specifies  Ts,  the  stopping  tolerance  for  sum  of  squares  con- 
vergence, i.e.,  for  convergence  based  on  the  relative  change  in  S{P,6).  The 
default  value  is  Ts  = where  ^ is  defined  as  the  smallest  value  such  that 
1 -f  ^ > 1 for  a <real>  computation  on  the  machine  being  used. 


SSTOL 

Ts 

< 0 
> 0 

- ^1/2 

= min{SST0L,  1} 

=>24  - PARTOL  [<real>  PARTOL] 

The  variable  that  specifies  Tp,  the  stopping  tolerance  for  parameter  conver- 
gence, i.e.,  for  convergence  based  on  relative  change  in  the  estimated  param- 
eters P and  A.  When  the  model  is  explicit  the  default  value  is  Tp  = and 
when  the  model  is  implicit  the  default  value  is  Tp  = where  ^ is  defined 
as  the  smallest  value  such  that  1 + ^ > 1 for  a <rGal>  computation  on  the 
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machine  being  used. 


PARTOL 

Tp 

< 0 
> 0 

= default  value 
= nim{PARTOL,  1} 

=>25  - MAXIT  [INTEGER  MAXIT] 

The  variable  that  specifies  7/,  the  maximum  number  of  iterations  allowed.  The 
default  value  depends  on  whether  the  fit  is  a restart  or  not.  (See  argument 
20-JOB.)  If  the  fit  is  not  a restart,  then 


T/  = 50  . 

If  the  fit  is  a restart,  then 

7/  = 7}_  + 10  , 

where  7/_  is  the  number  of  iterations  completed  in  the  previous  run,  thus 
indicating  that  the  procedure  will  continue  for  an  additional  10  iterations. 


MAXIT  Restart 

Ti 

<0  no 

yes 

>0  no 

yes 

= 50 

= Tj.  + 10 
= MAXIT 
= Tj-  + MAXIT 

If  MAXIT  = 0 then  no  iterations  wiU  be  taken,  but  whatever  computations  are 
required  to  complete  the  final  computation  report  will  be  made.  For  example, 
by  setting  MAXIT  = 0 and  the  third  digit  of  JOB  to  zero,  the  user  can  compute 
the  covariance  matrix  Vp  for  the  input  values  /3  and  A.  (See  arguments  20-JOB 
and  26-IPRINT,  and  also  §1.D.) 

=>26  - IPRINT  [INTEGER  IPRINT] 

The  variable  that  controls  generation  of  the  computation  reports  described  in 
§1.D.  The  default  computation  reports  include 

• a long  initial  summary, 

• no  iteration  summary,  and 

• a short  final  summary. 

When  IPRINT  < 0,  the  default  reports  are  generated  only  on  Ccr,  the  logical 
unit  specified  by  argument  28-LUNRPT.  When  IPRINT  > 0,  it  is  assumed 
to  be  a 4 digit  INTEGER  with  decimal  expansion  where  each  digit 

controls  a different  part  of  the  computation  report  and  whether  that  report 
is  to  be  generated  only  on  Ccr  or  to  both  Ccr  and  unit  6.  (See  §1.D). 
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Initial  summary 

14  - unit  - 

Ccr  6 

Iteration 
I3  Ja 

summary 
- unit  - 
Ccr  6 

Final  summary 

X\  - unit  - 

Ccr  6 

> 0 

= 0 

none 

none 

= 0 

none 

none 

= 0 

> 1 

none 

none 

= 0 

none 

none 

= 1 

short 

none 

= 1 

> 1 

short 

none 

= 1 

short 

none 

= 2 

long 

none 

= 2 

> 1 

long 

none 

= 2 

long 

none 

= 3 

short 

short 

= 3 

> 1 

short 

short 

= 3 

short 

short 

= 4 

long 

short 

= 4 

> 1 

long 

short 

= 4 

long 

short 

= 5 

short 

long 

= 5 

> 1 

short 

long 

= 5 

short 

long 

= 6 

long 

long 

= 6 

> 1 

long 

long 

= 6 

long 

long 

If  J2  = 0 no  iteration  summary  will  be  generated,  even  if  the  value  of  J3 
is  nonzero. 

If  J2  > 1 an  iteration  summary  will  be  generated  every  X2th  iteration  be- 
ginning with  iteration  one. 

=>27  - LUNERR  [INTEGER  LUNERR] 

The  variable  that  specifies  Cer,  the  logical  unit  number  to  be  used  for  error 
messages.  (See  §1.D.)  By  default,  Cer  = 


LUITERR 

Cer 

< 0 

= 6 

= 0 

error  messages  suppressed 

> 0 

= LUHERR 

=>28  - LUNRPT  [INTEGER  LUNRPT] 

The  variable  that  specifies  Ccr,  the  logical  unit  number  to  be  used  for  com- 
putation reports.  (See  also  argument  26-IPRINT,  and  §1.D.)  By  default, 
Ccr  = 6. 


LUHRPT 

Ccr 

< 0 

= 6 

= 0 

computation  reports  suppressed. 

even  when  argument  26-IPRIHT  is  nonzero 

> 0 

= LUHRPT 

=>29  - STPB  [<real>  STPB(NP)] 

The  singly  subscripted  array  that  specifies  the  relative  step  sizes,  K = 
1, . . .,p,  used  to  compute  the  finite  difference  derivatives  with  respect  to  each 
of  the  parameters  fS  as  discussed  in  §4.A.  The  default  value  is  set  as  described 
in  §4.A.iii  depending  on  whether  forward  or  central  finite  difference  derivatives 
are  being  computed.  (See  argument  20-JOB.) 
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STPB(l) 

II 

o o 

VI  A 

= default  value 
= STPB(K) 

=i^30  - STPD  [<real>  STPD(LDSTPD,M)] 

The  doubly  subscripted  array  that  specifies  the  relative  step  sizes,  /iau,  I = 
and  J = used  to  compute  the  finite  difference  derivatives 

with  respect  to  the  errors  in  each  of  the  elements  of  X as  discussed  in  §4.A. 
The  default  value  is  set  as  described  in  §4.A.iii  depending  on  whether  forward 
or  central  finite  difference  derivatives  are  being  computed.  (See  argument 
20-JOB.)  STPD  is  a structured  argument:  only  the  specific  elements  of  STPD 
identified  in  the  table  below  are  referenced  by  ODRPACK.  (See  §1.G.) 


STPDCl.l) 

LDSTPD 

For  I = 1, . . ., n, 

& J = 1, . . .,  m, 
hAu 

< 0 

— 

= default  value 

> 0 

= 1 

= STPD(1,J) 

> n 

= STPDd.J) 

=i^31  - LDSTPD  [INTEGER  LDSTPD] 

The  leading  dimension  of  array  STPD.  LDSTPD  must  either  exactly  equal  one, 
or  must  equal  or  exceed  n.  See  argument  30-STPD  for  further  details. 

=>32  - SCLB  [<real>  SCLB(NP)] 

The  singly  subscripted  array  that  specifies  the  scale  values,  SCALe{)3k}i 
K = 1, . . .,p,  of  the  function  parameters,  i.e.,  the  reciprocals  of  the  expected 
magnitudes  or  typical  sizes  of  — l,...,p.  Scaling  is  used  within  the 

regression  procedure  in  order  that  the  units  of  the  variable  space  will  have 
approximately  the  same  magnitude.  This  increases  the  stability  of  the  pro- 
cedure, but  does  not  affect  the  problem  specification.  Scaling  should  not  be 
confused  with  the  weighting  matrices  and  ws^  specified  by  arguments  11- 
WE  and  14-WD.  (See  §1.A,  and  §1.F.)  Except  as  noted  below,  the  scale  values 
specified  for  each  value  of  must  be  greater  than  zero.  The  default  values 
are  set  as  described  in  §4.D.i. 


SCLB(l) 

ForK=  l,...,p, 

scalb{/3k} 

o o 

VI  A 

= default  value 
= SCLB(K) 
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==>33  - SOLD  [<real>  SCLD(LDSCLD,M)] 

The  doubly  subscripted  array  that  specifies  the  scale  values  for  the  errors  A in 
the  explajiatory  variable  X,  i.e.,  the  reciprocals  of  the  expected  magnitudes  or 
typical  sizes  of  Au,  I = 1, . . . , n,  and  J = 1, . . . , m.  Scaling  is  used  within  the 
regression  procedure  in  order  to  ensure  that  the  units  of  the  variable  space  will 
have  approximately  the  same  magnitude.  This  increases  the  stability  of  the 
procedure,  but  does  not  affect  the  problem  specification.  Scaling  should  not 
be  confused  with  the  weighting  matrices  tyq  and  tuis^  specified  by  arguments 
11-WE  and  14-WD.  (See  §1.A,  and  §1.F.)  Except  as  noted  below,  the  scale 
values  specified  for  each  value  of  A must  be  greater  than  zero.  The  default 
values  are  set  as  described  in  §4.D.ii.  SOLD  is  a structured  argument:  only 
the  specific  elements  of  SOLD  identified  in  the  table  below  are  referenced  by 
ODRPACK.  (See  §1.G.) 


SCLDd.l) 

LDSCLD 

For  I = 1, . . .,  n, 

& J = 1, ..  .,m, 

scale{Au} 

< 0 

— 

= default  value 

> 0 

= 1 

= SCLDd.J) 

> n 

= SCLDd.J) 

=>34  - LDSCLD  [INTEGER  LDSCLD] 

The  leading  dimension  of  array  SOLD.  LDSCLD  must  either  exactly  equal  one, 
or  must  equal  or  exceed  n.  See  argument  33-SCLD  for  further  details. 

^35  - WORK  [<real>  WORK(LWORK)] 

The  singly  subscripted  array  used  for  <rGal>  work  spcice,  and  an  array  in 
which  various  computed  values  are  returned.  The  smallest  acceptable  dimen- 
sion of  WORK  is  given  below  in  the  definition  of  argument  36-LWORK.  The  work 
area  does  not  need  to  be  initialized  by  the  user  unless  the  user  wishes  to  ini- 
tialize A,  which  is  stored  in  the  first  n x m locations  of  WORK.  An  easy  way  to 
access  these  values,  either  for  initialization,  as  is  necessary  when  the  fourth 
digit  of  argument  20-JOB  is  nonzero  and  the  fit  is  by  orthogonal  distance  re- 
gression, or  for  analysis  upon  return  from  ODRPACK,  is  to  include  in  the 
user’s  program  the  declaration  statements 
<real>  DELTA «N>.<M» 

EQUIVALENCE  (WORK(l),  DELTA (1,1)) 

where  <N>  indicates  that  the  first  dimension  of  the  array  DELTA  must  be  ex- 
actly the  number  of  observations,  N = n;  and  <M>  indicates  that  the  second 
dimension  of  the  array  DELTA  must  be  exactly  the  number  of  columns,  M = m. 
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of  the  explanatory  variable  X.  This  allows  the  error  associated  with  observa- 
tion X(I,  J)  to  be  accessed  as  DELTA(I,J)  rather  than  as  WORKCl-t-CJ-D^N). 
The  values  in  DELTA  will  be  over  written  by  the  final  estimates  of  the  errors 
in  the  explanatory  variable  when  this  equivalencing  method  is  used.  Other 
values  returned  in  array  WORK  may  also  be  of  interest  and  can  be  accessed  as 
described  in  §5. A. 

N.B.,  if  the  fit  is  a “restart,”  i.e.,  if  the  fifth  digit  of  argument  20-JOB  is 
nonzero,  then  all  elements  of  vector  WORK,  including  the  values  of  DELTA,  must 
be  exactly  as  returned  from  a previous  call  to  ODRPACK. 

=>36  - LWORK  [INTEGER  LWORK] 

The  length  of  array  WORK. 

For  orthogonal  distance  regression  LWORK  must  equal  or  exceed 

18-|-llp-|-p^-fm-|-m^-|-4n9-l-6nm-|-2n9p-(-2n9m-|-5^-|-59-|-g(p-l-m)-f-(LDWE*LD2WE)g. 
For  ordinary  least  squares  LWORK  must  equal  or  exceed 

18  -I-  lip  -|-  p^  -I-  m -f  -f  Anq  4-  2nm  -f  2nqp  + 5^  4-  q{p  4-  m)  4-  (LDWE*LD2WE)9. 

^37  - IWORK  [INTEGER  IWORK(LIWORK)] 

The  singly  subscripted  array  used  for  INTEGER  work  space,  and  an  array  in 
which  various  computed  values  are  returned.  The  smallest  acceptable  dimen- 
sion of  IWORK  is  given  below  in  the  definition  of  argument  38-LIWORK. 

Certain  values  returned  in  array  IWORK  are  of  general  interest  and  can  be 
accessed  as  described  below  in  §5.B.  In  particular,  the  results  of  the  derivative 
checking  procedure  are  encoded  in  IWORK.  These  results  may  be  especially 
useful  if  ODRPACK’s  error  reports  have  been  suppressed.  (See  argument 
27~LUNERR.) 

N.B.,  if  the  fit  is  a “restMt,”  i.e.,  if  the  fifth  digit  of  argument  20-JOB  is 
nonzero,  then  all  elements  of  vector  IWORK  must  be  exactly  as  returned  from 
a previous  call  to  ODRPACK. 

=>38  - LIWORK  [INTEGER  LIWORK] 

The  length  of  array  IWORK. 

For  both  orthogonal  distance  regression  and  ordinary  least  squares  LIWORK 
must  equal  or  exceed 
20  -f-  p 4-  g(p  4-  m). 

<=39  - INFO  [INTEGER  INFO] 

The  variable  used  to  indicate  why  the  computations  stopped. 
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INFO 

Stopping  Condition 

= 1 

sum  of  squares  convergence 

= 2 

parameter  convergence 

= 3 

both  sum  of  squares  and  parameter  convergence 

= 4 

iteration  limit  reached 

> 5 

questionable  results  or  fatal  errors  detected 

When  INFO  > 5 the  questionable  results  or  fatal  errors  detected  by  ODR- 
PACK are  reported  in  the  messages  generated  on  the  logical  units  specified  by 
arguments  27-LUNERR  and  28-LUNRPT.  In  this  case,  INFO  is  a 5 digit  INTEGER 
with  decimal  expansion  where  J5  = 0 indicates  that  questionable 

conditions  were  found,  and  X5  > 1 indicates  that  fatal  errors  were  detected. 
The  nonzero  values  of  I5,  J4,  J3,  J2  and  X\  are  used  to  identify  what  condi- 
tions were  detected  at  the  time  the  program  stopped. 


Questionable  Results: 

J5  = 0 with  J4  ^ 0 
X3/O 

X2/O 
Ji  = l 
= 2 

= 3 

= 4 

derivatives  possibly  not  correct  (see  §4.A) 

ISTOP  0 at  last  call  to  FCN 

(see  argument  1-FCN) 

problem  is  not  full  rank  at  solution 

sum  of  squares  convergence 

parameter  convergence 

sum  of  squares  and  parameter  convergence 

iteration  limit  reached 

Falal  Errors: 

J5  = 1 with  J4  ^ 0 

N < 1 

M < 1 

NP  < 1 or  NP  > N 

JlT^O 

NQ  < 1 

J5  = 2 with  J4  ^ 0 

LDY  and/or  LDX  incorrect 

LDWE,  LD2WE,  LDWD  and/or  LD2WD  incorrect 

X2#0 

LDIFX,  LDSTPD  and/or  LDSCLD  incorrect 

Xi/0 

LWORK  and/or  LIWORK  too  small 

J5  = 3 with  J4  / 0 

STPB  and/or  STPD  incorrect 

SCLB  and/or  SOLD  incorrect 

X2  5^0 

WE  incorrect 

XiT^O 

WD  incorrect 

X5  = 4 

error  in  derivatives 

X5  = 5 

ISTOP  / 0 at  last  call  to  FCN 

II 

(see  argument  1-FCN) 
numerical  error  detected 
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Note  that  INFO  = 60000  indicates  an  error  possibly  caused  by  incorrectly 
specified  user  input  to  ODRPACK,  and  more  commonly  by  a poor  choice  of 
scale  or  weights,  or  a discontinuity  in  the  derivatives. 

2.C.  Examples 

The  following  sample  programs  invoke  DODR  and  DODRC  to  solve  the  examples  of  explicit, 
implicit  and  mtdtiresponse  problems  shown  in  §1.A. 

The  first  program  invokes  DODRC  with  user  supplied  derivatives,  the  second  program 
invokes  DODR  with  the  derivatives  approximated  by  ODRPACK  using  forward  finite 
differences,  and  the  third  invokes  DODRC  with  central  finite  difference  derivatives.  The 
use  of  forward  or  central  difference  derivatives  generally  causes  very  little  change  in  the 
results  from  those  obtained  using  analytic  derivatives.  (See  §4.A.) 

Users  are  encouraged  to  use  these  examples,  modified  as  necessary,  to  form  their  own 
ODRPACK  drivers.  Single  precision  sample  programs  can  be  easily  generated  from  these 
two  programs  by  changing  all  DOUBLE  PRECISION  variables  to  REAL,  and  substituting 
SODRfor  DODR  and  SODRC  for  DODRC.  Note  especially  that  by  using  NAXN,  MAXN,  HAXNP  and 
NAXNQ  to  specify  the  largest  problem  the  program  can  solve  without  modification,  and 
by  specifying  LWORK  and  LI  WORK  exactly  as  shown,  the  user  greatly  reduces  the  number 
of  changes  that  must  be  made  to  the  program  in  order  to  solve  a larger  problem. 


2.C.i.  Example  Problem  for  an  Explicit  Model 

The  following  sample  program  invokes  DODRC  to  solve  example  3.2.2  on  pages  230-238 
of  [Fuller,  1987].  The  data  (xi,yi)  3se  the  percent  saturation  of  nitrogen  gas  in  a brine 
solution  forced  into  the  pores  of  sandstone,  emd  the  observed  compressional  wave  velocity 
of  ultrasonic  signals  propagated  through  the  sandstone,  respectively.  These  data  (listed 
in  §2.C.i.b)  are  modeled  by  the  explicit  function 

Vi  « /(xi;/3)  = +)32[c^*‘  - l]^  i = l,...,n. 

The  starting  values  for  the  model  parameters  are 

P = (1500.0,-50.0,-0.1)'' 


and  A is  initialized  to  zero. 

Fuller  notes  that  it  is  reasonable  to  believe  that  the  saturation  measurements  of  0% 
and  100%  are  more  precise  than  the  other  saturation  measurements.  We  have  thus 
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“fixed”  xi,  X2  and  X12  at  their  original  values.  As  a consequence,  = ^|  = = 0 

at  the  solution.  We  assume  the  remaining  observed  data  axe  all  of  equal  precision,  and 
thus  set  Ttfe  = = 1 using  ODRPACK ’s  structured  argument  feature.  The  remaining 

arguments  are  set  to  their  default  values.  (See  §1.G.) 


2.C.i.a.  User  Supplied  Code 


PROGRAH  sanple 


ODRPACK 

<*=> 

<-»> 


Argvment 

Icn 


np 

nq 

beta 

y 

Idy 

z 

Idz 

we 

Idwe 

ld2«e 

«d 

Idvd 

ld2sd 

if  izb 

if  izz 

Idifz 

job 

ndigit 

tauf ac 

satol 

partol 

aazit 

iprint 

lunerr 

lunrpt 

atpb 

atpd 

Idstpd 

sclb 

acid 

Idscld 

vork 

Ivork 

iwork 


Definitions 

nane  of  the  user  si;q)plied  function  subroutine 
number  of  observations 

columns  of  data  in  the  ezplanatory  variable 

number  of  parameters 

number  of  responses  per  observation 

function  parameters 

response  variable 

leading  dimension  of  array  y 

ezplanatory  variable 

leading  dimension  of  array  z 

"epsilon"  weights 

leading  dimension  of  array  ve 

second  dimension  of  array  we 

"delta"  weights 

leading  dimension  of  array  wd 

second  dimension  of  array  wd 

indicators  for  "firing"  parameters  (beta) 

indicators  for  "firing"  ezplanatory  variable  (z) 

leading  dimension  of  array  if izz 

task  to  be  performed 

good  digits  in  subroutine  function  results 
trust  region  initialization  factor 
sun  of  sqiiares  convergence  criterion 
parameter  convergence  criterion 
mazimum  number  of  iterations 
print  control 

logical  unit  for  error  reports 

logical  unit  for  computation  reports 

step  sizes  for  finite  difference  derivatives  wrt  beta 

step  sizes  for  finite  difference  derivatives  wrt  delta 

leading  dimension  of  array  stpd 

scale  values  for  parameters  beta 

scale  values  for  errors  delta  in  ezplanatory  variable 

leading  dimension  of  array  scld 

DOUBLE  PRECISION  work  vector 

dimension  of  vector  work 

INTEGER  work  vector 
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c «>  livork  diBension  of  vector  ivork 

c <■«  info  stopping  condition 

c Paraneters  specifying  ■»T-inii  problen  sizes  handled  by  this  driver 
c Bam  Bari BUB  nuBber  of  observations 

c BeucB  BaziauB  ntuber  of  coluBns  in  explanatory  variable 

c Basip  ■wTiBiia  nuBber  of  function  paraneters 

c Baznq  BaziBUB  nuBber  of  responses  per  observation 

c Paraneter  Declarations  and  Specifications 

INTEGER  Idif z , Idscld , Idstpd , Idsd , Idve , Idz , Idy , ld2Bd , ld2Be , 

livork. Ivork, BazB,Bazn,Baziq>, Baznq 
PARAMETER  (Baza«5,Bazn’E25,BaznpBB,Baznq«l, 


BazB  *1’  BazBee2 


Idy^azn , ldz<>Bazn , 

Idsesl , ld2vesl , Idvd^l , ld2vd=l , 

Idif  z^azn , ldstpd=  1 , ldscld=l , 
lvork=18  llVBaznp  Baznpee2 

AeinTm^inTnq  -f  fiewmemaxa  + ?*iBitxn  ♦iBiiTnq»«iinmp  + 

2eBazn*BaznqVBazB  + Baznqee2  + 

SvBaxnq  Baznq*(Bazzq>-hnazB)  ldvevld2veVBaznq, 
livork=2(Haaz]:q>-^aznq*  (Baznp+Bazn)  ) 


Variable  Declarations 

INTEGER  i,info,iprint, j, job,l,lunerr,lunrpt,B,Bazit,n, 

+ ndigit,np,nq 

INTEGER  if izb(Baznp) ,if izz(ldifz,Bzum) , ivork (livork) 

DOUBLE  PRECISION  partol,sstol,tauf ac 

DOUBLE  PRECISION  beta(Baznp) ,sclb(Baznp) ,scld(ldscld,BazB) , 

+ stpb(Baznp) ,stpd(ldstpd,BazB) , 

+ vd (Idvd , ld2vd .Bazn) . ve ( Idve , ld2ve , Baznq) , 

vork(lvork),z(ldz,BazB)  ,y(ldy, Baznq) 

EXTERNAL  fen 


c Specify  default  values  for  dodre  argUBents 


ved.l.l) 

vdCl.l.l) 

ifizb(l) 

ifizzd.l) 

job 

ndigit 

taufac 

sstol 

partol 

Bazit 

iprint 

lunerr 

lunrpt 

stpbCl) 


-l.OdO 

-l.OdO 

-1 

-1 

-1 

-1 

-l.OdO 

-l.OdO 

-l.OdO 

-1 

-1 

-1 

-1 

-l.OdO 
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stpd(l.l)  - -l.OdO 
sclbCl)  - -l.OdO 
Bcldd.l)  “ -l.OdO 

c Set  up  ODRPACK  report  files 
lunerr  > 9 
liinrpt  « 9 

OPEH  (unit«9,f lies’ reportl*) 

c Read  problem  data,  and  set  nondefault  value  for  argument  ifixz 
OPEH  (unitsS.f iles>datal ’) 

READ  (5,FMTs*)  n,m,np,nq 
READ  (5,FMT-*)  (beta(i) ,i-l,np) 

DO  10  isi.n 

READ  (B.FMT-*)  (x(i, j),jsl,m) ,(y(i,l),l=l,nq) 
if  (x(i,l) .eq.O.OdO  .or.  x(i,l) .eq.lOO.OdO)  then 
ifixx(i,l)  s 0 

else 

ifixx(i,l)  s 1 
end  if 
10  CONTINUE 

c Specify  task:  explicit  orthogonal  distance  regression 
c «ith  user  supplied  derivatives  (checked) 

c covariance  matrix  constructed  with  recomputed  derivatives 

c delta  initialized  to  zero 

c not  a restart 

c and  indicate  short  initial  report 

c short  iteration  reports  every  iteration,  and 

c long  final  report 

job  s 00020 

iprint  > 1112 

c Compute  solution 
CALL  dodrcCfcn, 


+ 

n,m,np,nq. 

+ 

beta. 

+ 

y,ldy,x,ldx. 

+ 

se , Idse , ld2se , vd , Idvd , ld2vd  ^ 

+ 

if  ixb , if  ixx , Idif  X , 

+ 

j ob , ndigit , tauf  ac , 

+ 

sstol,partol,naxit , 

+ 

iprint , lunerr ,lunrpt , 

+ 

stpb, stpd, Idstpd, 

+ 

sclb , scld, Idscld, 

+ 

vork , Ivork , ivork , livork , 

+ 

END 

info) 
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STJBROUTIHE  lcn(n,B,np,nq, 

+ Idn.ldB.ldnp, 

+ bata,xplusd, 

+ ilixb.ilixx.ldifx, 

-4-  ideval,!  ,f  jacb,f  jacd, 

-t-  iatop) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


Subroutine  Arguments 

B=>  n number  of  observations 

n>  B number  of  columns  in  explanatory  variable 

»>  np  number  of  parameters 

nq  number  of  responses  per  observation 

■»>  Idn  leading  dimension  declarator  equal  or  exceeding  n 

»>  Idn  leading  dimension  declarator  equal  or  exceeding  a 

»>  Idnp  leading  dimension  declarator  equal  or  exceeding  np 

■b:>  beta  current  values  of  parameters 

»>  xplusd  current  value  of  explanatory  variable,  i.e.,  x delta 
ss>  ifizb  indicators  for  "fixing”  parameters  (beta) 

»=>  ifixx  indicators  for  "fixing"  explanatory  variable  (x) 

■°=>  Idifx  leading  dimension  of  array  ifixx 

»>  ideval  indicator  for  selecting  computation  to  be  performed 

<•**  f predicted  function  values 

<>«  f jacb  Jacobian  with  respect  to  beta 

<mm  f jacd  Jacobian  with  respect  to  errors  delta 
<"  istop  stopping  condition,  where 

0 means  current  beta  and  x+delta  were 
acceptable  and  values  were  computed  successfully 

1 means  current  beta  and  x+delta  are 

not  acceptable;  ODRPACK  should  select  values 
closer  to  most  recently  used  values  if  possible 
-1  means  current  beta  and  x+delta  are 
not  acceptable:  ODRPACK  should  stop 


Input  arguments,  not  to  be  changed  by  this  routine: 

INTEGER  i , ideval , ist  op , 1 , Idifx , Idm, Idn , Idnp ,a , n , np , nq 

DOUBLE  PRECISION  beta(np) ,xplusd(ldn,n) 

INTEGER  if  ixb(np)  ,ifixx(ldifx,n) 

Output  arguments: 

DOUBLE  PRECISION  f (ldn,nq) ,fjacb(ldn, Idnp, nq) ,fjacd(ldn, 1dm, nq) 
Local  variables 


INTRINSIC  exp 


c Check  for  imacceptable  values  for  this  problem 
IF  (beta(l)  .LT.  O.OdO)  THEN 
istop  = 1 
return 
ELSE 
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istop  E 0 
EHD  IF 


c Confute  predicted  values 

IF  (HOD(ideval.lO).GE.l)  THEN 
DO  110  1 « l.nq 
DO  100  i » l,n 

•=  beta(l)  + 

+ beta(2)*(exp(beta(3)*xplusd(i,l))  - 1.0d0)**2 

100  CONTINUE 

110  CONTINUE 
END  IF 


c 


Compute  derivatives  vith  respect  to  beta 
IF  (H0D(ideval/10.10).GE.l)  THEN 
DO  210  1 >=  l,nq 
DO  200  i > l.n 


f jacb(i,l,l) 
f jacb(i,2,l) 
ijacb(i,3,l) 

+ 

+ 

200  CONTINUE 

210  CONTINUE 
END  IF 


l.OdO 

(exp(beta(3)*xplusd(i,l))  - 1.0d0)**2 
beta(2)*2e 

(exp(beta(3)*xplusd(i,l))  - l.OdO)* 
exp (bet a (3 ) explusd ( i . 1) ) *xplusd ( i , 1 ) 


c Compute  derivatives  with  respect  to  delta 
IF  (M0D(ideval/100,10).GE.l)  THEN 
DO  310  1 l.nq 
DO  300  i B l.n 

fjacd(i,l,l)  K beta(2)*2* 

+ (ezp(beta(3)explusd(i,l))  - l.OdO)* 

+ exp(beta(3)*xplusd(i,l)}*beta(3) 

300  CONTINUE 

310  CONTINUE 
END  IF 


RETURN 

END 
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2.C.i.b.  User  Supplied  Data  (file  datal) 


12 

1 3 

1500.0 

-50.0  -( 

0.0 

1265.0 

0.0 

1263.6 

5.0 

1258.0 

7.0 

1254.0 

7.5 

1253.0 

10.0 

1249.8 

16.0 

1237.0 

26.0 

1218.0 

30.0 

1220.6 

34.0 

1213.8 

34.5 

1215.5 

100.0 

1212.0 
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2.C.i.c.  Report  Generated  by  ODRPACK  (file  reportl) 


*«*******«***«********«***************4i**«************* 

* odzpack  version  2.01  of  06-19-92  (double  precision)  * 

******************************************************* 


***  derivative  checking  report  for  fit  by  method  of  odr  *** 


' response 

1 of 

observation 

3 

user 

supplied 

relative 

derivative 

derivative 

wrt 

value 

difference 

assessment 

beta( 

1) 

l.OOD+00 

O.OOIHOO 

verified 

beta( 

2) 

1.S5D-01 

1.66D-06 

verified 

beta( 

3) 

1 . 19D+02 

2.94D-06 

verified 

delta ( 3. 

1) 

-2.39D+00 

2.96D-06 

verified 

number  of  reliable  digits  in  function  results  16 

(estimated  by  odrpack) 

number  of  digits  of  agreement  required  between 

user  supplied  and  finite  difference  derivative  for 

user  supplied  derivative  to  be  considered  verified  4 

row  number  at  which  derivatives  were  checked  3 

-values  of  the  explanatory  variables  at  this  row 

i(  3,  1)  B.OOOOOOOODfOO 

e«******«******«***e******4i****4i****e*e«eeeeee«e*4i**4i4i* 

* odrpack  version  2.01  of  06-19-92  (double  precision)  * 

******************************************************* 


***  initial  summary  for  fit  by  method  of  odr  *** 

problem  size: 


n 

nq 

m 

np 


12 

1 

1 

3 


(number  with  nonzero  weight 


(number  unfixed 


3) 


12) 


control  values: 

job  <B  00020 
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> abcda,  vhere 

a>0  »>  fit  is  not  a restart. 
b>0  »>  deltas  are  initialized  to  zero. 
c>0  »>  covariance  natriz  vill  be  computed  using 
derivatives  re-evaluated  at  the  solution. 
d^2  »>  derivatives  are  supplied  by  user, 
derivatives  sere  checked, 
results  appear  correct. 
e>0  »>  nethod  is  explicit  odr. 
ndigit  c 16  (estimated  by  odrpack) 

taufac  = l.OOD+00 


stopping  criteria: 

sstol  « 1.49D-08  (sum  of  squares  stopping  tolerance) 

partol  3.67D-11  (parameter  stopping  tolerance) 

mazit  ~ 50  (mazimum  number  of  iterations) 

initial  weighted  stub  of  squares  « 6.63720354D+05 

sum  of  squared  weighted  deltas  « 0 . OOOOOOOOD+OO 

sum  of  squared  weighted  epsilons  = 6 . 63720354D+05 

***  iteration  reports  for  fit  by  nethod  of  odr  *** 


it. 

num. 

cum. 

no.  fn 

evals 

weighted 

8um-of-sqs 

act.  ral. 
STUB-of-sqs 
reduction 

pred.  rel. 
sum-of-sqs 
reduction 

tau/pnorm 

g-n 

step 

1 

19 

2.51166IH01 

9.9996D-01 

9.9997D-01 

3.499D-01 

yes 

2 

20 

2.14730D+01 

1.4507D-01 

1.5249D-01 

6.900D-02 

yes 

3 

21 

2.144631H01 

1.2418D-03 

1.4631D-03 

8.573D-03 

yes 

4 

22 

2. 144551H01 

3.7676D-05 

4.7131D-05 

2.016D-03 

yes 

5 

23 

2. 14455D-i-01 

1.5116D-06 

1 . 8950D-06 

3.927D-04 

yes 

6 

24 

2.14455IH01 

6.2002D-08 

7.7723D-08 

8.017D-05 

yes 

7 

25 

2.14455D-t-01 

2.5359D-09 

3.1790D-09 

1.617D-05 

yes 

***  final  summary  for  fit  by  method  of  odr  *** 


— stopping  conditions: 

info  >■  1 »>  sum  of  squares  convergence, 

niter  > 7 (number  of  iterations) 

nfev  B 25  (number  of  function  evaluations) 

njev  > 8 (number  of  Jacobian  evaluations) 

irank  « 0 (rank  deficiency) 

rcond  > 6.12D-03  (inverse  condition  number) 

istop  > 0 (returned  by  user  from  subroutine  fen) 


final  weighted  sums  of  squares 


2.14455017IH01 
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8UB  of  squared  weighted  deltas  « 7.78974669IH00 
SUB  of  squared  weighted  epsilons  > 1.36557550D+01 


residual  standard  deviation 

degrees  of  freedom 


9 


1.54364294IH00 


estimated  beta(j) . j ■ 1,  ....  np: 


1 

beta 

1.2646B481IH-03 

s.d.  beta 

1.0349IH00 

9BX  confidence  interval  

1.26231139IH03  to  1 . 26699822D-K)3 

2 

-B.40184100IH01 

1.B840D-I-00 

-B.760B0942IH01 

to  -B.043172B7D+01 

3 

-8.78497122D-02 

6.3322D-03 

-1.02187862D-01 

to  -7.3B11B621D-02 

estimated  epsilonCi)  and 

deltaCi,*), 

i - 1,  . . . , n: 

i epsilonCi.l)  delta(i,l) 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 


-3.45194935D-01 
1 . 05480B06IH00 
-3.00719286D-02 
-1.13916405D-01 
-1.40250730D-01 
-5.53155556D-01 
-6.99564762D-01 
1.88412530IH00 
-1.70916306D+00 
1.80916198D+00 
1.90299896D-01 
-1.34707485D+00 


O.OOOOOOOOD+00 
0,00000000D+00 
-6.50838155D-02 
-2.67201445D-01 
-3.31357554D-01 
-1.30641313IH00 
-1.3252B687D+00 
1 .4B88B497D-)-00 
-1.18803B77D+00 
7.71243449D-01 
8.241392B3D-02 
0 . OOOOOOOOD+00 
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2.C.ii.  Example  Problem  for  an  Implicit  Model 

Tliis  sample  program  invokes  DODR  to  solve  the  implicit  problem  shown  in  example 
3.2.4  on  page  244  of  [Fuller,  1987].  In  this  example,  the  data  (listed  in  §2.C.ii.b)  are 
observations  digitized  from  the  x-ray  image  of  a hip  prosthesis,  where  the  variables 
Xi  = {vi,hi),  i = l,...,n,  are  the  vertical  and  horizontal  distances  from  the  origin, 
respectively,  and  the  implicit  model  is  that  of  the  ellipse 

fi{xi‘,/3)  = P3{vi-l3iy  + 2P^{vi-fii)ihi-/32)  + fi5{hi-fi2y  = 0 

for  T = 1, . . . ,n.  The  starting  values  for  the  model  parameters  are 

P = (-1.0, -3.0, 0.09, 0.02, 0.08) 


and  A is  initialized  to  zero.  Since  the  observed  data  are  all  of  equal  precision,  we  set 
vis  = 1 using  ODRPACK ’s  structured  argument  feature.  The  remaining  arguments  are 
set  to  their  default  values.  (See  §1.G.) 


2.C.ii.a.  User  Supplied  Code 


PROGRAM  sample 


ODRPACK  Argument 
=>  fen 
*=>  n 
»=>  m 
»=>  np 
=>  nq 
<«a=>  beta 

“>  y 

■=>  Idy 

•=>  I 

=>  Idx 
=>  we 
=>  Idwe 
■=>  ld2we 
=>  wd 
=>  Idwd 
=>  ld2wd 
=>  job 
■=>  iprint 
■*>  lunerr 
=>  lunrpt 
<==>  work 


Definitions 

name  of  the  user  supplied  function  subroutine 
number  of  observations 

columns  of  data  in  the  explanatory  variable 
number  of  parameters 
ntuber  of  responses  per  observation 
function  parameters 

response  variable  (unused  when  model  is  implicit) 

leading  dimension  of  array  y 

explanatory  variable 

leading  dimension  of  array  x 

initial  penalty  parameter  for  implicit  model 

leading  dimension  of  array  we 

second  dimension  of  array  we 

“delta"  weights 

leading  dimension  of  array  wd 

second  dimension  of  array  wd 

task  to  be  performed 

print  control 

logical  unit  for  error  reports 
logical  unit  for  computation  reports 
DOUBLE  PRECISION  work  vector 
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c 

c 

c 

c 


“>  Isork 
<“  ivork 
“>  1 ivork 
<“  info 


dinension  of  vector  vork 
INTEGER  vork  vector 
dimension  of  vector  ivork 
stopping  condition 


c 

c 

c 

c 

c 


Paraaeters  specifying  naxiauB  problea  sizes  handled  by  this  driver 
main  aaziaua  nuaber  of  observations 

aaza  Baziaua  nuaber  of  coluans  in  explanatory  variable 

aaxnp  maxiaua  nuaber  of  function  paraaeters 

iq  Baxiaua  nuaber  of  responses  per  observation 


Paraaeter  Declarations  and  Specifications 
INTEGER  Idvd,ldve,ldx,ldy,ld2vd,ld2ve, 

+ livork. Ivork, aaza, Baxn,Baxiq>,Baxnq 

PARAMETER  (BaxB=5,aaxn=25.aaznp=5,aaxnq>2, 
ldy=naxn,ldx>cBaxn, 

+ Idve=l,ld2ve<=l,ldvd=l,ld2vd=l, 

+ Ivork^lS  -f  lleaaxnp  Baxnpv*2  + aaxa  + BaxB**2  + 

+ 4eBaxneBaxnq  + SvBaxnVBaxB  + 2vBaxnVBajQiqVBaxnp  + 

+ 2eBaxnVBaxnqVBazB  + Banqee2  + 

■f  Svaaxnq  aaxnq* (aaxiy •^Baxa)  + ldve*ld2veVBajaiq, 

+ livork=20->-Baxz9+Baxnq*  (aaxi^+Baxa)  ) 


c Variable  Declarations 

INTEGER  i , inf  o , ipr  int . j , j ob , lunerr , lunrpt , a , n , np , nq 

INTEGER  ivork(livork) 

DOUBLE  PRECISION  betaCaaxnp), 

+ vd (Idvd , ld2vd .aaxa) , ve ( Idve , ld2ve .aaxnq) , 

-i-  vork(lvork)  ,x(ldx,BarB)  ,y(ldy, aaxnq) 

EXTERNAL  fen 


c 


Specify  default 
ve(l,l,l)  = 
vdd.l.l)  - 
job  •= 

iprint  = 
lunerr  >= 
lunrpt  = 


values  for  dodr  arguaents 
-l.OdO 
-l.OdO 
-1 
-1 
-1 
-1 


c Set  up  ODRPACK  report  files 
lunerr  > 9 

lunrpt  ~ 9 

OPEN  (unitB9,f ile=’report2') 


c Read  problea  data 

OPEN  (unites, file=*data2') 

READ  (5,FMT'>*)  n,B,np.nq 
READ  (B,FMT=*)  (beta(i) ,i=l,np) 
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DO  10  i-l.n 

READ  (B.FMT-*)  (i(i, j). j-l,*) 
10  CONTIHUE 


c 

c 

c 

c 

c 


Specify  task:  ii^licit  orthogonal  distance  regression 

with  forward  finite  difference  derivatiwes 
covariance  natrix  constructed  with  recomputed  derivatives 
delta  initialized  to  zero 
not  a restart 
job  ■ 00001 


c 


Compute 

solution 

CALL 

dodrCfcn, 

+ 

n,m,np,nq. 

+ 

beta. 

+ 

y,ldy,x,ldx. 

+ 

we , Idee , ld2we , wd , Idsd , ld2wd , 

+ 

job. 

+ 

ipr int , lunerr , lunrpt , 

+ 

work , Isork , iwork , liwork , 

+ 

END 

info) 

SUBROUTTNE  f cn(n,n,np,nq, 

+ ldn,ldm,ldnp, 

+ beta.zplusd, 

+ if ixb,if irx.ldifx, 

+ ideval.f ,f jacb.f jacd, 

+ istop) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Subroutine  Arguments 

a>  n number  of  observations 

mm>  ■ number  of  columns  in  explanatory  variable 

«>  np  number  of  parameters 

— > nq  number  of  responses  per  observation 

-=>  Idn  leading  dimension  declarator  equal  or  exceeding  n 

n>  1dm  leading  dimension  declarator  equal  or  exceeding  a 

=>  Idnp  leading  dimension  declarator  equal  or  exceeding  np 

n>  beta  current  values  of  parameters 

*=>  xplusd  current  value  of  explanatory  variable,  i.e.,  x delta 
a>  ifixb  indicators  for  "fixing”  parameters  (beta) 

ifixx  indicators  for  "fixing"  explanatory  variable  (x) 

=>  Idifx  leading  dimension  of  array  ifixx 

>=>  ideval  indicator  for  selecting  computation  to  be  performed 

<«<=  f predicted  function  values 

<=  f jacb  Jacobian  with  respect  to  beta 

f jacd  Jacobian  with  respect  to  errors  delta 
<s=  istop  stopping  condition,  where 
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c 0 Beans  current  beta  and  x+delta  sere 

c acceptable  and  values  were  computed  successfully 

c 1 Beans  current  beta  and  x+delta  are 

c not  acceptable:  ODRPACK  should  select  values 

c closer  to  Bost  recently  used  values  if  possible 

c -1  Beans  current  beta  and  x+delta  are 

c not  acceptable;  ODRPACK  should  stop 

c Input  arguBents,  not  to  be  changed  by  this  routine: 

IHTEGER  i.ideval,istop,l,ldifx,ldB,ldn,ldiq>,B,n,np,nq 

DOUBLE  PRECISION  beta(np) .xplusdCldn.m) 

INTEGER  ifixb(np) ,ifixx(ldifx,B) 

c Output  arguBents : 

DOUBLE  PRECISION  f (Idn.nq) .f jacb(ldn,l<inp.nq) ,f jacd(ldn.ldB.nq) 


c Check  for  unacceptable  values  for  this  problem 
IF  (beta(l)  ,GT.  O.OdO)  THEN 
istop  1 
return 
ELSE 

istop  > 0 
END  IF 


c 


Compute  predicted  values 

IF  (MOD(ideval.lO).GE.l)  THEN 
DO  110  1 = l.nq 
DO  100  i - l.n 

f (i,l)  K beta(3)*(xplusd(i,l)-beta(l))*e2  + 

+ 2ebeta(4)*(xplusd(i,l)-beta(l))* 

+ (xplusd(i,2)-beta(2))  + 

+ beta(5)*(xplusd(i,2)-beta(2))*e2  - l.OdO 


100  CONTINUE 

110  CONTINUE 
END  IF 


RETURN 

END 
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2.C.u.b.  User  Supplied  Data  (file  data2) 


20  2 

5 1 

-1.0  ■ 

-3.0  0.09  0.02  0.08 

0.50 

-0.12 

1.20 

-0.60 

1.60 

-1.00 

1.86 

-1.40 

2.12 

-2.54 

2.36 

-3.36 

2.44 

-4.00 

2.36 

-4.75 

2.06 

-5.25 

1.74 

-5.64 

1.34 

-5.97 

0.90 

-6.32 

-0.28 

-6.44 

-0.78 

-6.44 

-1.36 

-6.41 

-1.90 

-6.25 

-2.50 

-5.88 

-2.88 

-5.50 

-3.18 

-5.24 

-3.44 

-4.86 

54 


Using  ODRPACK 


2.C.ii.c.  Report  Generated  by  ODRPACK  (file  rGport2) 


****4t**********************************«*************** 

* odrpack  TerBion  2.01  oi  06-19-92  (double  precision)  * 


***  initial  Bumnary  lor  lit  by  nethod  ol  odr  *** 


problen  size: 


n “ 

20 

(number  with  nonzero 

weight 

nq  - 

1 

m > 

2 

np  - 

5 

(number  unlized  > 

5) 

control  values: 

job  = 00001 

« abode,  where 

a=0  «=> 
b-0  — > 
c=0  => 

d«0  -=> 
e=l  *=> 

ndigit  - 15 

taiilac  K 1 . OOD+OO 

stopping  criteria: 

sstol  ■:  1.49D-08  (stm  ol  squares  stopping  tolerance) 

partol  > 6.06D-06  (parameter  stopping  tolerance) 

mazit  > 100  (maziaun  number  ol  iterations) 


lit  is  not  a restart, 
deltas  are  initialized  to  zero, 
covariance  matriz  vill  be  computed  using 
derivatives  re-evaluated  at  the  solution, 
derivatives  are  estimated  by  lorvard  dillerences. 
method  is  implicit  odr. 

(estimated  by  odrpack) 


O.OOOOOOOOIHOO 

8.39823392D-01 
8.39823392D-01 
l.OD+01 

lunction  parameter  summary: 


initial  sum  ol  squared  weighted  deltas  ^ 

initial  penalty  lunction  value  >= 

penalty  term  > 

penalty  parameter  ■ 


indez  beta(k)  lized 


(k) 


(ilizb) 


scale  derivative 

step  size 
(sclb)  (stpb) 


1 -l.OOOOOOOOD+00 

2 -3.00000000D+00 

3 9.00000000D-02 

4 2.00000000I>-02 

6 8.00000000D-02 


no 

l.OOOOOOOOD+OO 

3.16228D-10 

no 

3.33333333D-01 

3.16228D-10 

no 

l.llllllllIHOl 

3.16228D-10 

no 

E.OOOOOOOOD+Ol 

3.16228D-10 

no 

1.25000000D+01 

3.16228D-10 
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explanatory  variable  and  delta  veight  auanary: 


index 

x(i.j) 

delta(i, j) 

fixed 

scale 

veight 

derivative 

(i.j) 

(if ixx) 

(scld) 

(wd) 

step  size 

(stpd) 

1.1 

5.000D-01 

O.OOOD+00 

no 

2.000+00 

1.000+00 

3.162280-10 

n,l 

-3.440D+00 

O.OOOD+00 

no 

2.910-01 

1 . 000+00 

3.162280-10 

1.2 

-1.200D-01 

O.OOOD+00 

no 

8.330+00 

1.000+00 

3.162280-10 

n.2 

-4.8600+00 

O.OOOD+00 

no 

2.060-01 

1.000+00 

3.162280-10 

***  final  aunnary  for  fit  by  nethod  of  odr  *** 


stopping  conditions: 

info  K 2 »>  paraneter  convergence, 
niter  > 18  (nunber  of  iterations) 

nfev  > 217  (niuber  of  function  evaluations) 

irank  « 0 (rank  deficiency) 

rcond  > 3.18D-02  (inverse  condition  number) 

istop  B 0 (returned  by  user  from  subroutine  fen) 


final  sum  of  squared  weighted  deltas  « 

final  penalty  function  value  = 
penalty  term  > 

penalty  parameter 


8 . 82420346D-02 

8. 824456 16D-02 
2 . 52700897D-06 
1 . OD+OS 


residual  standard  deviation 

degrees  of  freedom 


15 


7.66994283D-02 


estimated  beta(j) , j ■ 1,  ....  np: 


beta 


s.d.  beta 


95X  confidence  interval  


1 -9.99380972D-01 

2 -2.93104848D+00 

3 8.75730479D-02 

4 1 . 62299739D-02 

5 7.97538008D-02 


1.1138D-01 

1.0977D-01 

4.1061D-03 

2.7500D-03 

3.4963D-03 


-1.23682206IH00  to 
-3. 1650435 IIHOO  to 
7.88199915D-02  to 
1.03676338D-02  to 
7 . 23007073D-02  to 


-7.61939883D-01 

-2.69705344D+00 

9.63261044D-02 

2.20923140D-02 

8.72068944D-02 
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2.C.iii.  Example  Problem  for  an  Explicit  Model  with  Multiresponse  Data 


The  problem  shown  here  is  an  example  of  multiresponse  data  that  originates  because  the 
underlying  data  are  complex.  The  problem  is  described  in  Chapter  4,  and  on  pages  280- 
281,  of  [Bates  and  Watts,  1988].  In  this  case,  the  dependent  variable  is  the  pair  of  values 
representing  the  real  eind  imaginary  parts  of  complex-valued  impedance  measurements 
of  a polymer,  Zi,  i = 1, ...,n,  and  the  explanatory  variable,  Xj,  i = l,...,n,  is  the 
(real-valued)  frequency.  The  data  are  shown  in  §2.C.iii.b.  The  function  form  is  explicit, 
representing  the  dielectric  constant  by  the  general  model  proposed  by  [Havriliak  and 
Negami,  1967] 


for  i = l,...,n,  where  ] = 1.  For  ODRPACK,  this  must  be  encoded  as  two- 

term  multiresponse  data  with  yi  £ 8?^  representing  the  pair  of  values  [3?(2i),  5(2^)],  i = 
l,...,n.  Havriliak  and  Negami  (1967)  show  that  the  real  and  imaginary  components 
can  be  written  as 


where 


and 


9t{zi)  = /32  + {01  - 02)R^  cos{/3s<f>) 
= {01  - 02)R^  sm{Ps4>) 


R 


2 _ 


1 -I-  (27rxi/^3)^*  cos(7r/54/2)]  -f  [(2xi,/^3)^*  sin(x  M^)\ 
{2T:xil02tY*  sin(7r/34/2) 


(j)  = arctan 


1 + {2'KXil0zY^  cos(x/34/2) 


The  estimation  procedure  described  in  [Bates  ajid  Watts,  1988]  for  this  multiresponse 
problem  is  slightly  different  from  that  implemented  in  ODRPACK.  In  particular,  their 
procedure  provides  an  estimate  of  We,  but  does  not  include  estimates  of  A.  Thus,  we 
would  not  assume  that  the  results  obtained  here  using  ODRPACK  will  exactly  equal 
those  presented  by  Bates  and  Watts. 

For  our  example,  we  have  set  the  starting  values  for  the  model  parameters  to  be  the 
hnal  solution  shown  on  page  152  of  [Bates  and  Watts,  1988],  i.e.. 


0 = (4.398,2.451,8.245,0.487,0.571)’’ 


and  have  initialized  A to  the  decade  corrections  described  by  them. 

Bates  and  Watts  assume  that  there  is  no  error  in  the  first  decade  of  frequency  values, 
and  we  have  done  the  same,  “fixing”  these  variables  at  their  input  values.  Bates  and 


Using  ODRPACK 


57 


Watts  also  identify  two  outliers  in  the  data  set,  which  are  eliminated  from  our  analysis 
by  setting  the  corresponding  weights  to  zero.  The  remedning  weights  w^i  in  our  example 
are  set  to  an  estimate  of  the  2x2  covariance  matrix  of  the  errors  in  the  responses  of 
the  dependent  variable, 

where  Cj  is  the  estimate  of  i=  1, . . . ,n,  obtained  using  the  Bates  and  Watts  solution. 
The  weights  are  set  to  values  proportional  to  the  magnitude  of  the  frequencies.  The 
remaining  arguments  are  set  to  their  default  values.  (See  §1.G.) 


2.C.iii.a.  User  Supplied  Code 


PROGRAM  sanple 


c ODRPACK  Argiuent 

c “=>  Icn 

c ■*>  n 

C ■=>  B 

c “>  np 

c “>  nq 

c <**>  beta 

c -*>  y 

c =>  Idy 

c =>  X 

c ■=>  Idx 

c “>  ve 

c “>  Idwe 

c »>  ld2ve 

c ™>  vd 

c “>  Idwd 

c »>  ld2vd 

c “=>  ifixb 

c *=>  if  ixx 

c =>  Idifx 

c »=>  job 

c =>  ndigit 

c =«>  taufac 

c aatol 

c **>  partol 

c “>  Baxit 

c “>  iprint 

c B=>  lunerx 

c “>  lunrpt 

c =>  stpb 

c =>  stpd 

c =>  Idstpd 


Definitions 

nane  of  the  user  supplied  function  subroutine 
nuBber  of  observations 

coluBns  of  data  in  the  explanatory  variable 

nuBber  of  paraneters 

nuBber  of  responses  per  observation 

function  paraBeters 

response  variable 

leading  dinension  of  array  y 

explanatory  variable 

leading  diaension  of  array  x 

"epsilon"  weights 

leading  dinension  of  array  se 

second  dinension  of  array  we 

"delta"  weights 

leading  dinension  of  array  wd 

second  dinension  of  array  wd 

indicators  for  "fixing"  paraneters  (beta) 

indicators  for  "fixing"  explanatory  variable  (x) 

leading  dinension  of  array  ifixx 

task  to  be  perfomed 

good  digits  in  subroutine  fen  results 

trust  region  initialization  factor 

SUB  of  squares  convergence  criterion 

paraaeter  convergence  criterion 

BaxinuB  nunber  of  iterations 

print  control 

logical  unit  for  error  reports 

logical  unit  for  coaputation  reports 

step  sizes  for  finite  difference  derivatives  wrt  beta 

step  sizes  for  finite  difference  derivatives  wrt  delta 

leading  dinension  of  array  stpd 
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c 

c 

c 

c 

c 

c 

c 

c 


n>  8Clb 

“>  scld 
— > Idscld 
<“>  «ork 
“>  Ivork 
<“  ivork 
“>  1 ivork 
<*=  info 


scale  values  for  paraaeters  beta 

scale  values  for  errors  delta  in  explanatory  variable 

leading  dinension  of  array  scld 

DOUBLE  PRECISION  vork  vector 

dinension  of  vector  vork 

INTEGER  vork  vector 

dimension  of  vector  ivork 

stopping  condition 


c Parameters  specifying  maximum  problem  sizes  handled  by  this  driver 


c 

maxn 

maximum  number 

of 

observations 

c 

maxm 

maximum  number 

of 

columns  in  explanatory  variable 

c 

maxnp 

maximum  number 

of 

function  parameters 

c 

maxnq 

■aximun  number 

of 

responses  per  observation 

c 


Parameter  Declarations  and  Specifications 

INTEGER  Idif X , Idscld , Idstpd , Idvd , Idve , Idx , Idy , ld2vd , ld2ve , 

livork, Ivork, naxm,naxn,naxi^),naxnq 
PARAMETER  (maxm=5,meucn=100,maxiq>=25,maxnq=5. 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


ldy=maxn , Idx^maxn , 

ldve=maxn, ld2 ve=maxnq , ldvd=naxn , ld2vd=l , 
ldifx=maxn,ldscld=l,ldstpd=l , 

lvork=18  llenaxnp  + maxnpe*2  + naxa  aaza*e2  + 

4*naxn*naxnq  + 6*maxnemaxm  ^ennfm^imiTnqanimTTip  + 
2emaxnemaxnq*naxa  + maxnq*e2  + 

Eemaxnq  maznq*(aaxiip-fmaxm)  ldveeld2vevinaxnq, 
livork=2(H-maziq>4maxnq*  (maxnp-^maxm)  ) 


Variable  Declarations 

INTEGER  i,info,iprint,  j,  job,l,lunerr,lunrpt,n,nnxit,n, 

ndigit.np.nq 

INTEGER  if ixb(maxnp) ,ifizx(ldifx,maxm) , ivork (livork) 

DOUBLE  PRECISION  partol,sstol,tauf ac 

DOUBLE  PRECISION  beta(aaxnp) ,sclb(aaxnp) ,8cld(ldscld,maxm) , 
8tpb(naxnp),8tpd(ldstpd,maxm}, 

+ vd(ldvd,ld2vd,aaxn)  ,ve(ldve,ld2ve,maxnq) , 

+ vork(lvork) ,x(ldx,naxm) ,y(ldy,maxnq) 

EXTERNAL  fen 


c Specify  default  values  for  dodre  arguments 
ve(l,l,l)  = -l.OdO 
vdd.l.l)  = -l.OdO 
ifixb(l)  = -1 
ifixxd.l)  = -1 
job  = -1 

ndigit  = -1 
taufac  = -l.OdO 
sstol  = -l.OdO 
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partol 

m 

-l.OdO 

maxit 

m 

-1 

iprint 

m 

-1 

lunerr 

m 

-1 

lunrpt 

m 

-1 

stpb(l) 

m 

-l.OdO 

stpd(l,l) 

m 

-l.OdO 

sclb(l) 

m 

-l.OdO 

scld(l.l) 

B 

-l.OdO 

c Set  up  ODRPACK  report  files 
lunerr  « 9 

lunrpt  > 9 

OPEN  (uiiit=9,file®’report3’) 

c Reed  problem  data 

OPEN  (unit-5,lile-‘data3') 

READ  (5,FHT>*)  n,B,np,nq 
READ  (5,FMT=*)  (beta(i) ,i=l,np) 

DO  10  i=l,n 

READ  (S.FMT-*)  (x(i. j) , j=l.B) . (y(i,l) .I'l.nq) 

10  CONTINUE 

c Specify  task  as  explicit  orthogonal  distance  regression 
c with  central  difference  derivatives 

c covariance  matrix  constructed  with  recomputed  derivatives 

c delta  initialized  by  user 

c not  a restart 

c and  indicate  long  initial  report 
c no  iteration  reports 

c long  final  report 

job  - 01010 

iprint  » 2002 

c Initialize  delta,  and  specify  first  decade  of  frequencies  as  fixed 
DO  20  i^l.n 

if  (x(i,l) .It. lOO.OdO)  then 
vork(i)  B O.OdO 
ifixz(i,l)  ■■  0 

else  if  (x(i,l) .le.150.0d0)  then 
vork(i)  B O.OdO 
if ixx(i, 1)  B 1 

else  if  (x(i,l).le.l000.0d0)  then 
Bork(i)  B 25.0d0 
if ixx(i, 1)  B 1 

else  if  (x(i,l) .le. 10000. OdO)  then 
vorkCi)  B 560. OdO 
if ixx(i, 1)  B 1 

else  if  (x(i,l) .le. 100000. OdO)  then 
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■ork(i)  - 9500. OdO 
ifixxCi.l)  - 1 

else 

«ork(i)  - 144000. OdO 
ifizx(i,l)  - 1 
end  if 
20  CONTIVUE 


Set  Heights 
DO  30  i-l,n 
il 


He(i,l,l) 

He(i.l,2) 

se(i,2,l) 

ve(i,2.2) 

else 

se(i,l,l) 
se(i.l.2) 
ve(i.2.1) 
wed, 2, 2) 
end  il 
wd(i,l.l) 

30  COHTINOE 


(x(i,l) .eq. 100. OdO 
O.OdO 
O.OdO 
O.OdO 
O.OdO 


.or.  x(i,l) ■•q* IBO.OdO)  then 


559. 6d0 
-1634. OdO 
-1634. OdO 
8397. OdO 

(1.0d-4)/(x(i.l)**2) 


c Coiq>ute  solution 
CALL  dodrcClcn, 


+ 

n,n,np,nq. 

+ 

bets. 

+ 

7,ldy,x,ldx, 

+ 

we , Idwe , ld2we , vd, Idvd, ld2Hd, 

+ 

ilixb.if ixx.ldilx. 

+ 

j ob , ndigit , taul sc , 

+ 

sstol,partol,naxit , 

+ 

iprint , lunerr , Ixuirpt , 

+ 

stpb , stpd, Idstpd, 

+ 

sclb , sold, Idscld, 

+ 

work , Isork , isork , 1 isork , 

+ 

EHD 

info) 

SXXBROUTINE 

lcn(n,n,np,nq. 

+ 

ldn,ldm,ldnp. 

+ 

beta,xplusd. 

+ 

il ixb , il ixx , Idilx , 

+ 

ideval . 1 , 1 j acb , 1 j acd , 

+ 

istop) 

c Subroutine  Arguments 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


■=>  n 
— «>  ■ 

»=>  np 
m=>  nq 
«=>  Idn 
— > Ida 
««=>  Idnp 
*=>  beta 
*=>  xplusd 
=>  ilizb 
-=>  ifixi 
■*>  Idilx 
“*>  ideval 
<==  f 
<=  Ijacb 
<»=  Ijacd 
<==  iatop 


nuabar  ol  obserrationa 

nuaber  of  coluana  in  explanatory  variable 

nuaber  of  paraaetera 

nuaber  of  reaponsea  per  obaervation 

leading  diaenaion  declarator  equal  or  exceeding  n 

leading  diaenaion  declarator  equal  or  exceeding  a 

leading  diaenaion  declarator  equal  or  exceeding  np 

current  valuea  of  paraaetera 

currant  value  of  explanatory  variable,  i.e.,  x + delta 
indicatora  for  "fixing”  paraaetera  (beta) 
indicatora  for  "fixing"  explanatory  variable  (x) 
leading  diaenaion  of  array  if ixx 

indicator  for  aelecting  coaputation  to  be  perforaed 

predicted  function  values 

Jacobian  with  respect  to  beta 

Jacobian  with  respect  to  errors  delta 

stopping  condition,  where 

0 aeans  current  beta  and  x+delta  were 
acceptable  and  values  sere  coaputed  successfully 

1 aeans  c\irrent  beta  and  x+delta  are 

not  acceptable:  ODRPACK  should  select  values 
closer  to  aost  recently  used  values  if  possible 
-1  aeans  current  beta  and  x+delta  are 
not  acceptable;  ODRPACK  should  stop 


c Input  arguaents,  not  to  be  changed  by  this  routine: 

IlfTEGER  i , ideval , ist  op , Idif  x , Ida , Idn , Idnp  ,a , n , iq> , nq 

DOUBLE  PRECISION  beta(np) ,xplusd(ldn,n) 

INTEGER  if ixb(np) ,if ixx(ldifx,n) 

c Output  arguaents : 

DOUBLE  PRECISION  f (ldn,nq) ,f jacb(ldn,ldnp,nq) ,f jacd(ldn,ldB,nq) 
c Local  variables 

double  precision  freq, pi, onega,ctheta,stheta, theta, phi, r 
INTRINSIC  atan2,e:q>,sqrt 


c Check  for  unacceptable  values  for  this  problea 
DO  10  i>l,n 

IF  (xplusd(i,l).LT.0.0d0)  THEN 
iatop  * 1 
return 
END  IF 
10  CONTINUE 
iatop  0 

pi  3.141592653589793238462643383279d0 

theta  ~ pi*beta(4)*0.5d0 
ctheta  B cos (theta) 
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stheta  B 8in(theta) 

c Coiq>uta  predicted  values 

IF  (MOD(ideval.lO).GE.l)  THEN 
DO  100  1^1,11 

freq  > zplu8d(i,l) 

oaega  > (2.0d0«pi*lreq*exp(-beta(3)))**beta(4) 
phi  ■■  atan2((oBegae8theta} , (l+oaega*ctheta)) 
r “ (beta(l)-beta(2))  • 

+ 8qrt((l-i-oBega*ctheta)**2+ 

+ (oBega*8theta)**2)**(-beta(5)) 

f (i,l)  = beta(2)  + r*co8(beta(5)*phi) 
l(i,2)  = r*8in(beta(B)*phi) 

100  CONTINUE 
END  IF 

RETURN 

END 
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2.C.iii.b.  User  Supplied  Data  (file  dataS) 


23  1 

5 

2 

4.0  2.0  7 

.0  0.40 

0.50 

30.0 

4.220 

0.136 

50.0 

4.167 

0.167 

70.0 

4.132 

0.188 

100.0 

4.038 

0.212 

150.0 

4.019 

0.236 

' 200.0 

3.956 

0.257 

300.0 

3.884 

0.276 

500.0 

3.784 

0.297 

700.0 

3.713 

0.309 

1000.0 

3.633 

0.311 

1500.0 

3.540 

0.314 

2000.0 

3.433 

0.311 

3000.0 

3.358 

0.305 

5000.0 

3.258 

0.289 

7000.0 

3.193 

0.277 

10000.0 

3.128 

0.255 

15000.0 

3.059 

0.240 

20000.0 

2.984 

0.218 

30000.0 

2.934 

0.202 

50000.0 

2.876 

0.182 

70000.0 

2.838 

0.168 

100000.0 

2.798 

0.153 

150000.0 

2.759 

0.139 

64 


Using  ODRPACK 


2.C.iii.c.  Report  Generated  by  ODRPACK  (file  reports) 

» odrpack  Tersion  2.01  of  06-19-92  (double  precision)  * 

**^^f*4f************************************************* 


***  initial  sunmary  lor  lit  by  method  of  odr  *** 


problea  size: 


n ■ 

23 

(number  with  nonzero 

weight 

nq  - 

2 

m > 

1 

np  - 

5 

(number  unfized  > 

5) 

control  values: 

job  = 01010 

K abode,  where 

a«0  ■=> 
b=l  ■=> 
c»0  “> 

d-1  — > 

esO  -=> 

ndigit  15 

taufac  = I.OOIHOO 

stopping  criteria: 

sstol  B 1.49D-08  (s\in  of  squares  stopping  tolerance) 

partol  3.67D-11  (parameter  stopping  tolerance) 

mazit  50  (maziaum  number  of  iterations) 


fit  is  not  a restart, 
deltas  are  initialized  by  user, 
covariance  matriz  will  be  computed  using 
derivatives  re-evaluated  at  the  solution, 
derivatives  are  estimated  by  central  differences, 
method  is  ezplicit  odr. 

(estimated  by  odrpack) 


initial  weighted  sun  of  squares 

sum  of  squared  weighted  deltas 
sun  of  squared  weighted  epsilons 


1.71064070D+03 

2.01382943D-04 

1.71064050IH03 


function  parameter  summary: 


indez 

beta(k) 

f ized 

scale 

derivative 

(k) 

(if izb) 

(sclb) 

step  size 
(stpb) 

1 

4.00000000IH00 

no 

2.50000000D-01 

1.000000-05 

2 

2.000000001>f00 

no 

5 . OOOOOOOOD-01 

1.000000-05 

3 

7.00000000IH00 

no 

1.42857143D-01 

1.000000-05 

4 

4.00000000D-01 

no 

2 . 500000000+00 

1.000000-05 

5 

5.00000000D-01 

no 

2.000000000+00 

1.000000-05 

ezplanatory  variable  and  delta  weight  summary: 
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index 

x(i.j) 

delta(i, j) 

fixed 

scale 

veight 

derivative 

(i.j) 

(iiixx) 

(scld) 

(vd) 

step  size 
(stpd) 

1.1 

3.000IH01 

0.0000+00 

yes 

3.330-02 

1.110-07 

1.000000-05 

n,l 

1.500IHOB 

1.4400+05 

no 

6.670-06 

4.440-15 

1.000000-05 

response 

variable  and 

epsilon  error  veight 

summary: 

index 

(i.l) 

y(i.i) 

veight 

(ve) 

1.1 

4.2200400 

5 . 5960+02 

n,l 

2. 7590+00 

5 . 5960+02 

1.2 

1.3600-01 

8.3970+03 

n,2 

1.3900-01 

8.3970+03 

***  final  susBary  lor  lit  by  aathod  ol  odr  *** 


- — atopping  conditions: 

inlo  c 1 sun  ol  squares  convergence, 
niter  « 8 (number  ol  iterations) 

nlev  B 121  (number  ol  function  evaluations) 

irank  ^ 0 (rank  deficiency) 

rcond  ^ 8.1BD-03  (inverse  condition  number) 

istop  K 0 (returned  by  user  from  subroutine  Icn) 


final  veighted  sums  ol  squares  > 

sum  ol  squared  veighted  deltas  « 
sum  ol  squared  veighted  epsilons  = 


4.20538922D-01 

5.54021897D-04 

4.19984900D-01 


residual  standard  deviation 

degrees  ol  freedom 


16 


1.62122431D-01 


estimated  beta(j) . j « 1,  ....  np: 


beta 


s.d.  beta 


■ 95%  confidence  interval 


1 

4.379988030+00 

1.30630-02 

4.352293880+00 

to 

4.407682180+00 

2 

2.433305760+00 

1.30500-02 

2.405638200+00 

to 

2.460973320+00 

3 

8.002884590+00 

1.16710-01 

7.755448030+00 

to 

8.250321150+00 

4 

5.101147160-01 

1.32640-02 

4.819928240-01 

to 

5.382366090-01 

5 

5.173902330-01 

2.88530-02 

4.562184980-01 

to 

5.785619680-01 

estimated  epsilon(i)  and  delta(i,e),  i > 1, 


n: 


i epsilon(i,l)  epsilon(i,2) 


delta(i,l) 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
23 


-7.3855879BD-03 
-1.05614733D-03 
-2.70863920D-03 
4.68593517D-02 
8.08102389D-03 
1.53882522D-03 
4.60535703D-03 
4.50906164D-03 
-1.00621895D-03 
1.05810802D-02 
6.93622739D-03 
3. 9582801 lD-05 
-3.77617796D-03 
-5.56734978D-04 
2.08263807D-03 
-7.50689916D-03 
-1.56731844D-03 
-5.93223183D-04 
1 . 15260099D-04 
2. 6364111 lD-04 
-3.81011180D-04 
-3. 368226 llD-04 
2.87173883D-03 


1.26939187D-03 
-1.22846292D-03 
-2.14347329D-03 
-4.25940138D-03 
-3.47639194D-03 
3. 852937 13D-04 
1.19118896D-03 
1.235708920-03 
-2. 918650430-04 
3.272841940-03 
2.434821060-03 
1.759050140-05 
-2.429078140-03 
-1.701237840-03 
-2.237232330-03 
2 . 164628930-03 
2.033670850-04 
2.720691710-05 
-2.421261310-07 
5.185103190-06 
-1.039638500-05 
-1.261413910-05 
1.411998410-04 


0 . 000000000^00 
0.000000000-K)0 
0.000000000+00 
0.000000000+00 
0.000000000+00 
3.036944000+01 
3.789867500+01 
6.226304870+01 
1.111869800+02 
1 . 157098770+02 
2.414365910+02 
9.613445320+02 
1.330298450+03 
2.075115660+03 
2.902895320+03 
5.218158180+03 
7.545646360+03 
1.742010210+04 
2.427454720+04 
3.784920520+04 
5.534932800+04 
8.757914320+04 
1.294963000+05 


3.  When  the  Model  Is  Very  Time  Consuming 


ODRPACK  executes  user  supplied  subroutine  FCN  not  only  to  compute  the  initial  sum 
of  the  squared  errors  5(/3,  and  to  obtain  function  and  derivative  values  within  its  main 
iterative  procedure,  but  also  when  setting  the  default  value  for  the  number  of  good  digits 
in  the  function  results,  when  performing  derivative  checking,  and  when  constructing  the 
covariance  matrix  and  standard  deviations  of  the  estimated  parameters  When  the 
time  required  for  hnding  the  solution  is  dominated  by  the  evaluation  of  FCN,  the  user 
will  want  to  make  judicious  use  of  these  options  in  light  of  their  “cost.”  Let  p be  the 
number  of  unfixed  parameters  P,  and  let  y?  = 1 if  the  fit  is  by  orthogonal  distance 
regression  and  y?  = 0 if  the  fit  is  by  ordinary  least  squares.  Then  the  number  of  times 
the  function  and  derivatives  are  evaluated  in  each  of  these  instances  is  summcirized  as 


follows: 

Computation 

Function 

Derivative 

Controlling 

Evaluations 

Evaluations^ 

Variable 

• Initial  S{P,6): 

1 

0 

— 

• Per  Iteration: 

> 1 

1 

— 

• Default  number  of 

good  digits  in  function  results: 

4 

0 

NDIGIT 

• Derivative  checking: 

> q{p  + (pm) 

1 

JOB 

• Default  covariance  matrix: 

0 

1 

JOB 

Users  with  a very  time  consuming  subroutine  FCN  should  also  be  aware  of  two  of  ODR- 
PACK’s  options  that  are  specifically  designed  for  such  problems.  The  most  importeint 
of  these  is  the  restart  facility.  The  other  is  the  option  of  constructing  the  covariance 
matrix  without  recomputing  the  derivative  matrices  at  the  solution.  This  second  option 
is  discussed  in  §2.Bii,  under  the  description  of  subroutine  argument  JOB,  and  also  in 
§4.B.  The  remainder  of  this  section  describes  how  ODRPACK ’s  restart  facility  can  be 
used  to  minimize  the  risk  of  losing  important  results  because  system  imposed  time  limits 
are  reached  before  the  solution  is  found. 

forward  finite  difference  approximation  to  the  derivative  requires  9(p  + ^m)  function  evaluations, 
and  a central  finite  difference  approximation  requires  2q{p  + fm)  function  evaluations.  (See  §4.A.) 
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The  restart  facility  enables  the  user  to  step  through  the  solution  procedure  one  or  more 
iterations  at  a time  without  incurring  any  additional  function  or  derivative  evaluations 
over  what  would  be  required  if  the  procedure  were  allowed  to  run  to  convergence.  Figure 
3.1  shows  an  example  of  how  the  restart  facility  can  be  employed.  This  example  allows 
up  to  a total  of  30  iterations,  and  writes  the  contents  of  arrays  BETA,  WORK  and  IWORK 
to  a file  between  every  iteration,  alternating  between  two  files.  This  minimizes  the 
chance  of  losing  significant  amounts  of  important  information  due  to  system  imposed 
limits:  if  such  a limit  is  reached  before  convergence,  arrays  BETA,  WORK  and  IWORK  can 
be  restored  using  the  saved  data  and  the  computations  restarted.  In  this  example, 
the  initial  computation  report  is  only  generated  at  the  first  iteration,  while  the  final 
computation  report  is  generated  after  every  iteration.  The  options  selected  include 
constructing  the  covariance  matrix  using  the  derivative  matrices  from  the  last  iteration, 
and  thus  no  additional  calls  to  subroutine  FCN  are  incurred  in  order  to  provide  the 
standard  deviations  of  the  parameters  printed  in  each  of  the  reports. 

Users  with  very  time  consuming  problems  should  be  aware  that,  depending  on  which 
options  are  selected,  ODRPACK  will  make  at  least  one  and  possibly  more  calls  to 
subroutine  FCN  before  attempting  to  generate  any  computation  reports.  Also,  on  many 
systems  the  output  generated  by  a program  is  not  written  directly  in  a file  but  rather 
is  stored  in  a “buffer”  until  the  buffer  is  full,  and  is  only  written  to  the  file  at  that 
point.  If  a run  is  aborted  prematurely,  either  by  the  user  or  because  a system  imposed 
limit  is  reached,  then  the  content  of  these  buffers  might  not  be  emptied  into  the  files 
associated  with  them.  Thus,  the  files  associated  with  the  logical  units  specified  by 
arguments  LUNERR  and  LUNRPT  might  not  include  all  information  actually  generated  by 
ODRPACK  at  the  time  it  stopped.  When  the  user  is  not  getting  the  expected  reports 
from  ODRPACK,  it  may  be  necessary  to  have  ODRPACK  generate  reports  directly  to 
“standard  output,”  which  is  usually  not  buffered,  in  order  to  determine  exactly  where 
the  computations  are  stopping.  (See  §2.B.ii,  subroutine  arguments  IPRINT,  LUNERR  and 
LUNRPT.) 
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Figure  3.1:  Using  ODRPACK’s  Restart  Facility 


I 

c set  up  riles  to  ssts  cosqnitatioiis  for  fntnre  restarts 
lanl  B 11 

OPEI  (onitBlnnl ,f ile^ ’savel .dat ’ ) 
lim2  ■ 12 

OPEI  (uiitBlaii2,file<°’saTe2.dat’) 

c set  the  maximi  nraber  of  iterations  for  each  call  to  ODBPACK  to  one 
c so  results  can  be  stored  beteeen  iterations 
Mzit  ■ 1 

c set  argnaent  appropriately,  Mking  sure  for  first  iteration  that 
c fit  is  not  a restart  and 

c covariance  natriz  is  constructed  sithont  reconpating  derivatives 

job  * 00100 

c set  iprint  to  indicate  a long  initial  report, 
c a short  iteration  report,  and 

c a long  final  report 

iprint  - 2112 

c step  throngh  up  to  30  iterations 
DO  100  niter  « 1 , 30 

CALL  dodrc(fcn,  ...  ,info) 

c save  the  contents  of  beta,  eork  and  isork  for  future  reference 
IF  (■od(niter,2) .eq.l)  THEl 
lun  lunl 
ELSE 

lun  B lnn2 
EID  IF 

OPEI  (UIIT-lun) 

WRITE  (Inn,*)  (beta(k) ,kBl,np) 

WRITE  (Inn,*)  (Bork(i),isi,luork) 

WRITE  (lun,*)  (ivork(i) ,iBl,liBork) 

CLOSE  (UlIT-lun) 

IF  ( inf o.ge. 10000  .or.  ■od(info,10) .le.3)  then 
c stop  because  either  a fatal  error  uas  detected,  or  the  problem  converged 
stop 
ELSE 

c set  job  to  indicate  the  next  iteration  is  a restart  and 
c set  iprint  to  suppress  the  initial  report  for  future  iterations, 
job  B 10100 
iprint  * 0112 
EID  IF 
100  COITIIUE 
c 
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4.  Computational  Details 


4. A.  Computing  the  Jacobian  Matrices 


As  was  noted  in  §1.A,  the  matrices  of  first  partial  derivatives,  i.e.,  the  Jacobian  matrices 


k=l p,  &/  = l 

j = &/=!,. ..,9, 


(4.1) 


are  required  at  every  iteration.  These  can  be  provided  by  the  user  as  described  for 
subroutine  argument  FCN  in  §2.B.ii,  or  can  be  approximated  automatically  by  ODR- 
PACK.  User  supplied  derivatives  are  generally  either  “hand  coded”  as  is  done  in  the 
example  program  shown  in  §2.C.i,  or  are  the  product  of  an  “automatic  differentiation” 
tool.  ODRPACK’s  approximations  cire  formed  using  either  forward  or  central  finite 
differences. 


4.A.i.  “Hand  Coded”  Derivatives 

Hand  coded  derivatives  are  those  produced  by  the  user  without  the  aid  of  a differenti- 
ation tool.  Because  coding  errors  axe  a common  problem  with  hand  coded  derivatives, 
ODRPACK  has  an  option  to  check  the  validity  of  the  user  supplied  derivative  code  by 
comparing  its  results  to  finite  difference  values  for  the  derivative.  The  derivative  check- 
ing procedure  examines  the  unfixed  variables  at  only  one  row  of  the  Jacobian  matrix, 
and  is  therefore  quite  efficient.  Checking  only  one  row  is  reasonable  for  regression  mod- 
els since  the  same  code  is  frequently  used  to  compute  the  model  function  and  derivatives 
for  each  row,  as  is  the  case  for  each  of  the  examples  shown  in  §2.C. 

When  the  value  of  the  user  supplied  derivative  disagrees  with  the  corresponding  finite 
difference  value,  the  checking  procedure  attempts  to  determine  whether  the  disagree- 
ment is  due  to  an  error  in  the  user’s  code,  or  is  due  to  the  inaccuracy  of  the  finite 
difference  approximation.  The  checking  procedure  generates  an  error  report  when  one 
or  more  of  the  derivatives  are  found  to  be  questionable.  This  information  is  also  returned 
to  the  user  in  subroutine  argument  IWORK.  (See  §5.B.) 
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Questionable  derivatives  can  occur  when  the  derivative  is  exactly  zero  or  when  the  nu- 
merical derivative  used  in  the  checking  procedure  is  believed  to  be  inaccurate  because 
of  the  properties  of  the  function.  Zero  valued  derivatives  are  questionable  because  they 
could  indicate  that  the  initial  values  of  the  function  parameters  might  be  hiding  an 
error  in  the  derivative,  such  as  could  occur  if  the  initial  value  of  one  of  the  parameters 
was  zero.  Users  should  examine  the  ODRPACK  error  reports,  or  the  encoded  values  in 
subroutine  argument  IWORK  as  described  in  §5.B,  to  determine  the  cause  of  the  ques- 
tionable results,  and  then  examine  subroutine  FCN  to  insure  that  there  is  not  an  error  in 
the  user  supplied  derivatives  that  could  be  adversely  affecting  the  least  squares  results. 


4.A.ii.  Automatic  Differentiation 

Automatic  differentiation  tools  produce  code  to  calculate  a function’s  derivatives  di- 
rectly from  the  code  used  to  compute  the  function  values.  Such  tools  enable  the  user  to 
generate  the  derivatives  required  by  ODRPACK  without  the  tedium  and  errors  associ- 
ated with  hand  coded  derivatives.  An  overview  of  automatic  differentiation  is  presented 
in  [Griewank,  1989],  and  a survey  of  automatic  differentiation  software  is  provided  in 
[Juedes,  1991].  We  have  found  tools  such  as  ADIFOR  [Bischof  et  ai,  1991]  and  DAPRE 
[Stephens  and  Pryce,  1991],  which  are  precompilers  that  transform  a Fortran  subroutine 
that  evaluates  the  function  into  a Fortran  subroutine  that  evaluates  both  the  function 
and  its  derivatives,  especially  suitable  for  use  with  ODRPACK. 

Currently,  most  differentiation  tools,  including  ADIFOR  and  DAPRE,  generate  code  to 
evaluate  the  function  and  derivative  values  simultaneously.  Least  squares  procedures, 
however,  need  the  derivatives  only  after  determining  that  a satisfactory  new  point  has 
been  found,  and  this  determination  requires  at  least  one  and  possibly  more  function 
evaluations.  Thus,  when  using  a differentiation  tool,  one  has  two  choices:  either  the 
function  and  derivatives  can  be  always  evaluated  together  using  the  code  generated  by 
the  differentiation  tool;  or  the  hand  coded  function  can  be  evaluated  by  itself  until  a 
satisfactory  new  point  has  been  found,  at  which  time  the  derivative  code  generated 
by  the  differentiation  tool  code  can  be  evaluated.  The  first  choice  has  the  drawback 
that  sometimes  the  computed  derivative  values  will  not  be  used.  If  the  second  option 
is  selected,  the  function  values  produced  by  the  differentiation  tool  generated  code  are 
used  only  for  the  evaluation  of  the  Jacobian  matrices  and  not  within  the  least  squares 
procedure. 

For  the  differentiation  tools  and  problems  examined  by  the  authors,  the  time  required 
to  evaluate  the  derivative  and  function  together  using  code  generated  by  automatic 
differentiation  is  frequently  significantly  more  than  that  required  to  evaluate  only  the 
hand  coded  function.  We  thus  believe  that  in  most  cases  it  will  be  more  cost  effective 
for  users  to  employ  the  second  of  the  two  choices  mentioned  above.  Users  must  also 
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be  aware  that  ODRPACK  never  asks  that  both  the  function  and  its  derivatives  be 
computed  in  a single  call  to  subroutine  FCN,  and  that  it  is  possible  that  an  invocation 
of  FCN  for  evaluating  the  derivative  will  not  be  immediately  preceded  with  a call  to  FCN 
to  evaluate  the  function  for  the  same  parameter  values.  Thus,  if  the  first  of  the  two 
options  is  employed,  the  user  will  need  to  construct  a mechanism  for  saving  the  computed 
derivatives  until  ODRPACK  actually  requests  them,  possibly  using  the  Fortran  COMMON 
facility,  and  also  a mechanism  for  determining  whether  the  saved  derivative  values  were 
in  fact  those  evaluated  at  the  selected  parameter  values. 


4.A.iii.  Finite  Difference  Derivatives 

Finite  difference  derivatives  are  automatically  constructed  by  ODRPACK  to  approx- 
imate the  JacobicLn  matrices  when  the  user  does  not  supply  code  within  subroutine 
FCN  to  compute  them.  Either  forward  or  central  finite  differences  can  be  employed  for 
the  approximation.  (See  §2.B.ii,  subroutine  argument  JOB.)  A central  finite  difference 
derivative  gives  a more  accurate  approximation  than  the  corresponding  forward  finite 
difference  derivative,  but  at  the  expense  of  an  additional  call  to  subroutine  FCN  for  each 
partial  derivative  computed.  The  interested  reader  is  referred  to  [Dennis  and  Schnabel, 
1983]  and  [Gill  et  aL,  1981]  for  a more  complete  discussion  of  forward  and  central  finite 
difference  approximations. 


4.A.iii.a.  Forward  Finite  Difference  Derivatives 

The  forward  finite  difference  derivative  with  respect  to  /3  for  response  I of  observation  t 
is  computed  using 

fujxj  -H  Sj;  0 + h0^Uk)-  fujxj  + 6i;P) 

hu 

i = 1,. . .,n,  k = 1, . . .,p,  and  / = 1, . . .,g,  where  u*  is  the  kth.  unit  vector,  i.e.,  the  kth 
column  of  a p X p identity  matrix,  eind  where  is  the  finite  difference  step  size, 

hk  = A:  = 1, . . .,p,  (4.2) 

with  the  relative  step  size  for  parameter  0k  specified  by  subroutine  argument  STPB. 
(See  §2.B.ii.)  The  default  value  for  the  relative  step  for  a forward  finite  difference 
derivative  is 

where  V’  indicates  the  number  of  good  digits  in  the  results  of  the  user  supplied  subroutine 
FCN.  (See  §2.B.ii,  subroutine  argument  NDIGIT.)  This  default  value  is  selected  based  on 
empirical  evidence  that  indicates  it  generally  outperforms  the  commonly  recommended 
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value  h = (See,  e.g.,  [Dennis  and  Schnabel,  1983].)  Procedures  for  selecting 

near  optimal  relative  step  sizes  are  discussed  in  [Gill  et  al,  1981]  and  [Schnabel,  1982].) 

The  step  must  be  large  enough  so  that  approximately  half  of  the  good  di^ts  of 
faixi  + P + h^^Uk)  and  fa{xi  + Si; P)  will  be  the  same.  The  forward  finite  difference 
approximation  to  the  derivative  can  then  be  expected  to  have  roughly  half  the  number 
of  good  digits  as  are  in  the  computed  value  of  fa{xi  + Si;^).  When  the  computation  of 
/t/(®»  + ^»»/5)  has  sufficient  precision,  then  forward  finite  difference  derivatives  will  cause 
very  little  change  in  the  results  from  those  that  would  be  obtained  using  hand  coded  or 
automatic  differentiation  derivatives. 

The  forward  finite  difference  derivatives  with  respect  to  A are  formed  analogously  to 
those  with  respect  to  (3.  (See  §2.B.ii,  subroutine  argument  STPD,  for  specification  of  the 
relative  stepsize.) 

4.A.iii.b.  Central  Finite  Difference  Derivatives 

When  the  user  suspects  that  the  forward  finite  difference  approximation  will  not  provide 
sufficient  precision,  then  a central  finite  difference  approximation  can  be  used.  (See 
§2.B.ii,  subroutine  argument  JOB.)  The  central  finite  difference  approximation  to  the 
partial  derivative  with  respect  to  /S*  for  response  I of  observation  t is  given  by 

fujxi  + 6i;P  + hkUk)  - fiijxj  -\-Si;(3-  hkUk) 

hk 

t = 1, . . .,n,  fc  = 1, . . .,p,  and  I = The  step  hk  is  formed  as  in  (4.2),  where  the 

default  value  for  the  relative  step  for  a central  finite  difference  derivative  is 

h = 10"^/^  . 

The  central  finite  difference  derivatives  with  respect  to  A are  again  formed  analogously. 
4.B.  Covariance  Matrix 

The  linearized  confidence  regions  and  intervals  for  the  unknowns  and  A estimated 
by  orthogonal  distance  regression  are  the  same  as  the  linearized  regions  and  intervals 
that  would  be  obtained  if  the  orthogonal  distance  regression  problem  were  solved  as  a 
p+nm  parameter  nonlinear  ordinary  least  squares  problem.  If  we  express  the  orthogonal 
distance  regression  problem  defined  by  (1.8)  or  (1.9)  as  such  a nonlinear  ordinary  least 
squares  problem  with  nq  + nm  observations  and  p + nm  unknowns,  and  we  designate 


Computational  Details 


75 


the  unknowns  of  this  ordinary  least  squares  problem  as  rj'^  = (/?’',  SJ, . . . , 6^),  then  the 
sum  of  squares  to  be  minimized  is 

S{v)  = G{vr!lG(r,) 

where  G{ti)  is  the  vector  valued  function  whose  ith  “element”  is  defined  by 

/ \ / /t(®»  0)  ~ Vi  i = 1, . . . , Tl, 

i = n+l 2n. 

and  n 6 ^{rui+nm)x{nq+nm)  jg  block  diagonal  weighting  matrix  given  by 


n 


^1 


We. 


The  ordinary  least  squares  representation  of  (1.8)  or  (1.9)  is  thus 


2n 


min5(»7)  = min  5^^.(77)’^n.i(7,(7/) 

V V 

t=l 


(4.3) 


where  flu  denotes  the  (t,i)th  “element”  of  fl. 

Let  G'{Tj)  6 denote  the  Jacobian  matrix  with  {u,k)th.  element  equal 

to  dg„{T})fdTik  evaluated  at  K we  assume  that  G'{r})  and  ft  are  full  rank,  so  that 
[G'{‘q)'^flG'{‘q)]  is  nonsingiilar,  then  the  linearized  covariance  matrix  for  the  estimators 
fj  is  the  (p  + nm)  X (p  + nm)  matrix 

V = a^[G\vyflG'm-\ 

where  a = S{‘q)ffi  = S{0,6)lfi  is  the  estimated  residual  variance  with  p degrees  of 
freedom.  (The  degrees  of  freedom  is  the  number  of  observations  with  nonzero  weights 
minus  the  number  of  parameters  actually  being  estimated,  i.e.,  fi  = n — p.)  This 
covariance  matrix  V can  be  partitioned 


V = 


% %6  1 

V60  Vs  J 


where  G is  the  covariance  matrix  for  the  estimators  0,  Vs  E is  the 

covariance  matrix  for  the  estimators  A,  and  V^s  = Vg^  6 gives  covariances  be- 

tween 0 and  A.  It  is  the  covariance  matrix  V^  of  the  estimators  0 that  is  automatically 
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provided  by  ODRPACK.  The  actual  computational  technique  used  by  ODRPACK  to 
compute  Vp  is  described  in  detail  in  [Boggs  and  Rogers,  1990b]. 

By  default,  ODRPACK  will  recompute  the  Jacobian  matrices  at  the  Anal  solution  before 
constructing  the  covariance  matrix.  However,  ODRPACK  also  provides  the  option  of 
constructing  the  covariance  matrix  using  the  Jacobian  matrices  from  the  last  iteration. 
(See  §2.B.ii,  subroutine  argument  JOB.)  The  option  of  using  the  Jacobian  matrices 
from  the  last  iteration  to  construct  the  covariance  matrix  is  especially  useful  when  the 
evaluation  of  user  supplied  subroutine  FCN  is  very  time  consuming.  Assuming  that  the 
algorithm  has  actually  converged,  using  the  Jacobian  matrices  from  the  last  iteration 
should  give  essentially  the  same  covariance  matrix  as  that  which  would  be  obtained 
using  the  Jacobian  matrices  recomputed  at  the  solution.  Once  the  user  conArms  that 
the  solution  is  satisfactory,  the  covariance  matrix  can  easily  be  computed  at  the  actual 
solution  by  calling  ODRPACK  with  the  Anal  values  of  (3  and  A as  input,  and  with 
argument  MAXIT  = 0 and  the  third  digit  of  argument  JOB  set  to  0. 

The  steindard  deviations,  dp,  of  the  function  parameters  listed  in  the  ODRPACK  Anal 
report  are  the  square  roots  of  the  diagonal  elements  of  Vp,  i.e., 

% = V^^\k,k)  . 

The  95%  conAdence  intervals  are  computed  using 

Pk  ± <.975,/i^/3* 

where  <.975,^  is  the  appropriate  value  for  constructing  a two-sided  95%  conAdence  interval 
using  the  Student’s  t value  for  fj,  degrees  of  freedom.  When  p.  > 20,  t.975,^  « 2;  when 
p < 5,  <.975, /I  > 2.5  . 

K necessary,  the  full  covariance  matrix  V for  all  of  the  estimators  rj  can  be  computed 
using  the  equations  given  in  [Boggs  and  Rogers,  1990b],  or  can  be  “automatically” 
obtained  from  most  ordinary  least  squares  software  (including  ODRPACK)  by  solving 
the  orthogonal  distance  regression  problem  as  the  ordinary  least  squares  problem  deAned 
by  (4.3). 

Note  that  for  nonlinear  ordinary  least  squares,  the  linearized  conAdence  regions  and 
intervals  are  asymptotically  correct  as  n — > 00  [Jennrich,  1969].  For  the  orthogonal 
distance  regression  problem,  they  have  been  shown  to  be  asymptotically  correct  as  c*  — ♦ 
0 [Fuller,  1987].  The  difference  between  the  conditions  of  asymptotic  correctness  can  be 
explained  by  the  fact  that,  as  the  number  of  observations  increases  in  the  orthogonal 
distance  regression  problem  one  does  not  obtain  additional  information  for  A.  Note  also 
that  V is  dependent  upon  the  weight  matrix  D,  which  must  be  assumed  to  be  correct, 
and  cannot  be  con  Armed  from  the  orthogonal  distance  regression  results.  Errors  in  the 
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values  of  and  ujy,.  that  form  Cl  will  have  an  adverse  affect  on  the  accuracy  of  V and 
its  component  parts.  The  results  of  a Monte  Carlo  experiment  examining  the  accuracy 
of  the  linearized  confidence  intervals  for  four  different  measurement  error  models  is 
presented  in  [Boggs  and  Rogers,  1990b].  Those  results  indicate  that  the  confidence 
regions  and  intervals  for  A are  not  as  accurate  as  those  for  p. 

Despite  its  potential  inaccuracy,  the  covariance  matrix  is  frequently  used  to  construct 
confidence  regions  and  intervals  for  both  nonlinear  ordinary  least  squares  and  mea- 
surement error  models  because  the  resulting  regions  and  intervals  are  inexpensive  to 
compute,  often  adequate,  and  familiar  to  practitioners.  Caution  must  be  exercised  when 
using  such  regions  and  intervals,  however,  since  the  validity  of  the  approximation  will 
depend  on  the  nonlinearity  of  the  model,  the  variance  and  distribution  of  the  errors, 
and  the  data  itself.  When  more  reliable  intervals  and  regions  are  required,  other  more 
accurate  methods  should  be  used.  (See,  e.g.,  [Bates  and  Watts,  1988],  [Donaldson  and 
Schnabel,  1987],  and  [Efron,  1985].) 


4.C.  Condition  Number 


For  a linear  least  squares  system  of  equations 


AP=2Y  , 


(4.4) 


with  A € assumed  to  have  full  column  rank  and  =2  meaning  “equals  in  the  least 

squares  sense,”  the  condition  number  of  A is  defined  as 


K 


(A)  ^ ||44|1 


A^ 


where  A^  = {A'^A)~^A'^  is  known  as  the  pseudo  inverse  of  A.  From  this  definition,  we 
can  show  that  «(>!)  > 1,  and  that  k{A)  — ♦ 00  as  the  columns  of  A become  dependent. 

Using  /c(>l),  bounds  can  be  constructed  on  the  relative  error  in  the  true  least  squares 
solution  P*  = A^Y  due  to  a perturbation  Ey  of  Y , or  to  a perturbation  Ea  of  A. 
While  the  actual  bounds,  discussed  in  detail  in  [Stewart,  1973],  are  quite  complicated, 
we  can  roughly  approximate  them  as  follows.  Let  py  and  Pa  denote  the  solutions  to 
the  perturbed  systems,  py  = A^{Y  + Ey)  and  Pa  = {A  + Ea)^Y,  respectively,  and  let 
R = AP*  — Y denote  the  residual  at  P*.  Then 


\\P*-Py\\  < 

m\  "" 


\\p* 


(A) 

^ llrll 

(4.5) 

■HA)  H'l 

11^4111 

(4.6) 

and 
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For  (4.5)  we  thus  observe  that  the  relative  error  in  the  solution  could  be  as  much  as  /c(i4) 
larger  than  the  relative  error  in  Y.  Similarly,  for  (4.6)  we  observe  that  if  the  residual 
R is  small,  then  the  relative  error  in  A is  multiplied  by  ^(^4),  whereas  if  the  residual  is 
not  small,  then  the  second  term  in  (4.6)  will  dominate  and  the  relative  error  in  A could 
be  multiplied  by  as  much  as  k^{A). 

If  we  express  the  condition  number  k{A)  as  a power  of  10,  i.e.,  k{A)  = 10^,  then  (4.5) 
implies  that  the  elements  of  the  least  squares  solution  could  have  K.  fewer  significant 
digits  of  accuracy  than  the  elements  of  Y",  while  (4.6)  implies  that  the  least  squares 
solution  could  have  as  many  as  2K,  fewer  significant  digits  of  accuracy  than  the  elements 
of  A.  Therefore,  k{A)  sufficiently  large  could  indicate  that  the  least  squares  solution  has 
no  significant  digits.  For  the  condition  number  to  provide  a meaningful  estimate  of  ill- 
conditioning,  however,  the  user’s  problem  must  be  formulated  so  that  the  errors  in  the 
columns  of  A are  equilibrated.  This  requires  an  intimate  knowledge  of  the  problem,  and 
cannot  be  done  as  part  of  ODRPACK’s  automatic  scaling  procedure.  If  this  equilibration 
has  not  been  done  by  the  user,  the  condition  number  returned  by  ODRPACK  may  not 
reflect  the  true  conditioning  of  the  problem. 

ODRPACK  returns  an  approximation  to  the  inverse  of  the  condition  number  k(M^  J^), 
where  is  the  Jacobian  matrix  of  partial  derivatives  with  respect  to  fi,  and  Ms  is  a 
block  diagonal  matrix  formed  using  the  partial  derivatives  with  respect  to  A as  described 
in  [Boggs  et  al.,  1987].  The  matrix  MsJp  is  used  by  ODRPACK  to  form  a linearization 
of  the  user’s  problem  at  the  solution,  and  can  thus  be  substituted  for  the  matrix  A in 
in  the  above  discussions.  The  approximate  inverse  condition  number  is  calculated  as 
described  in  the  Linpack  Users’  Guide  [Dongarra  et  al.,  1980]. 

4.D.  Scaling  Algorithms 

Poorly  scaled  problems,  i.e.,  problems  in  which  the  unknowns  fi  and  A vary  over  several 
orders  of  magnitude,  can  cause  difficulty  for  least  squares  procedures.  ODRPACK’s 
scaling  algorithms  (discussed  below)  attempt  to  overcome  these  difficulties  automati- 
cally, although  it  is  preferable  for  the  user  to  choose  the  units  of  fi  and  A so  that  their 
estimated  values  will  have  roughly  the  same  magnitude.  (See,  e.g.,  [Dennis  and  Schn- 
abel, 1983].)  When  the  variables  have  roughly  the  same  magnitude,  the  ODRPACK 
scaling  algorithm  will  select  scale  values  that  are  roughly  equal,  and  the  resulting  com- 
putations will  be  the  same  (except  for  the  effect  of  finite  precision  arithmetic)  as  an 
unsealed  analysis,  i.e.,  an  analysis  in  which  all  of  the  scale  values  are  set  to  one.  If  the 
user  does  not  do  this,  the  ODRPACK  scaling  algorithm  will  select  varying  scale  values. 
This  will  not  change  the  optimal  solution,  but  it  may  affect  the  number  of  iterations 
required,  or,  in  some  cases,  whether  the  algorithm  is  or  is  not  successful. 
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The  scale  value  times  the  corresponding  absolute  value  of  the  expected  solution  should 
be  approximately  one.  For  example,  if  /3k  is  expected  to  be  of  order  10'°  then  scale{/3k} 
should  be  set  to  10”'°,  while  if  /3k  is  expected  to  lie  between  —10”^  and  —10”^  then 
scale{/3k}  should  be  set  to  10^.  Scaling  should  not  be  confused  with  the  weighting 
matrices  and  wsi  specified  by  subroutine  arguments  WE  and  WD.  (See  also  §1,A  cind 

§1-F.) 
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4.D.i.  Scaling  P 

ODRPACK  chooses  the  default  scale  values  for  the  estimated  values  of  as  follows. 

If  some  of  the  values  of  0 are  nonzero  then 

Let  Prnax  = the  largest  nonzero  absolute  value  in  0,  and 
Let  /Jniin  = the  smallest  nonzero  absolute  value  in  0. 

For  K = 1, ...  ,p  do 

If  = 0 then 

SCALE{/3k}  = lOfPmm 
Else  if  log(/3xnax)  - log(/3inin)  > 1 then 

scale{/3k}  = 1/\0k\ 

Else 

SCALE{/3k}  = l/0nuix 
Else  if  all  of  the  values  of  0 are  zero  then 
For  K = 1, . . .,p  do 

scale{;3k}  = 1 

Users  may  also  substitute  their  own  scaling  values  for  0.  (See  §2.B.ii,  subroutine  argu- 
ment SCLB.) 
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4.D.ii.  Scaling  A 

ODRPACK  chooses  default  scale  values  for  the  estimated  errors  Au  in  the  explanatory 
variables  as  follows. 

For  J = 1, ... ,m  do 

If  some  values  in  column  J of  X are  nonzero  then 

Let  Xniax  = the  largest  nonzero  absolute  value  in  column  J 
Let  Xnun  = the  smallest  nonzero  absolute  value  in  column  J 
For  I = l,...,n  do 

If  Xu  = 0 then 

SCALE{Au}  = lO/Xxnin 
Else  if  log(XnM«)  - log(Xnun)  > 1 then 
scale{Au}  = l/|Xu| 

Else 

scale{Au}  = 

Else  if  all  values  in  column  J of  X are  zero  then 
For  I = 1, . . .,n  do 

scale{Au}  = 1 

Users  may  also  substitute  their  own  scaling  values  for  A.  (See  §2.B.ii,  subroutine  argu- 
ment SOLD.) 
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5.  Work  Vectors 


5. A.  Extracting  Information  from  Vector  WORK 

Upon  return  from  a call  to  ODRPACK,  array  WORK  contains  values  that  may  be  of  inter- 
est to  the  user.  To  extract  information  from  WORK,  the  following  declaration  statement 
must  be  added  to  the  user’s  program: 

LOGICAL 
+ ISODR 
INTEGER 

+ DELTAI , EPSI , XPLUSI , FNI , SDI , VCVI , 

+ RVARI , WSSI , WSSDEI .WSSEPI .RCONDI ,ETAI , 

+ OLMAVI , TAUI , ALPHAI , ACTRSI , PNORMI , RNORSI , PRERSI , 

+ PARTLI , SSTOLI , TAUFCI , EPSMAI , 

+ BETAOI , BETACI , BETASI , BETANI . SI , SSI , SSFI , QRAUXI , UI , 

+ FSI,FJACBI,WE1I,DIFFI, 

+ DELTSI ,DELTNI ,TI .ITI , OMEGAI .FJACDI , 

+ WRKII . WRK2I , WRK3I , WRK4I , WRK5I , WRK6I , WRK7I , 

+ LWKMN 

where  DELTAI  through  WRK7I  are  variables  that  indicate  the  starting  locations  within 
WORK  of  the  stored  values,  and  LWKMN  is  the  minimum  acceptable  length  of  array  WORK. 

The  appropriate  values  of  DELTAI  through  HRK7I  are  obtained  by  invoking  subroutine 
SWINF  when  using  either  of  the  single  precision  ODRPACK  subroutines  SODR  or  SODRC, 
and  by  invoking  DWINF  when  using  either  of  the  double  precision  subroutines  DODR  or 
DODRC.  The  call  statements  for  SWINF  and  DWINF  have  the  same  argument  lists.  To 
invoke  either  subroutine,  use 

CALL  <winf> 

+ (N,M,NP,NQ,LDWE,LD2WE, ISODR, 

+ DELTAI , EPSI , XPLUSI , FNI , SDI , VCVI , 

+ RVARI , WSSI , WSSDEI , WSSEPI , RCONDI ,ETAI , 
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+ OLMAVI , TAUI , ALPHAI , ACTRSI , PNORMI , RNORSI , PRERSI , 

+ PARTLI , SSTOLI , TAUFCI , EPSMAI , 

+ BETAOI , BETACI , BETASI , BETANI , SI . SSI , SSFI , QRAUXI ,UI . 

+ FSI,FJACBI,WE1I.DIFFI, 

+ DELTSI .DELTNI ,TI ,TTI , OMEGAI , F JACDI , 

+ WRKII , WRK2I , WRK3I , WRK4I , WRK5I , WRK6I , WRK7I , 

+ LWKMN) 

where  SWINF  should  be  substituted  for  <winf  > when  using  SODR  and  SODRC,  and  DWINF 
should  be  substituted  for  <winf  > when  using  DODR  and  DODRC.  The  variables  N,  M NP, 
NQ,  LDWE,  and  LD2WE  must  be  input  to  SIWINF  and  DIWINF  with  exactly  the  same  values 
as  were  used  in  the  original  call  to  ODRPACK,  and  ISODR  must  be  input  set  to  true  if 
the  ht  was  by  orthogonal  distance  regression,  and  be  input  set  to  false  if  the  fit  was  by 
ordinary  least  squares.  Note  that  when  ISODR  is  false,  the  locations  that  are  specified  by 
DELTASI  through  WRKII  are  the  same  as  the  location  specified  by  DELTAI. 

In  the  following  descriptions  of  the  information  returned  in  WORK,  > indicates  values 
that  are  most  likely  to  be  of  interest. 

> WORK  (DELTAI)  is  the  first  element  of  an  n x m array  DELTA  containing  the  estimated 

errors  A in  the  explanatory  variables  at  the  solution,  where 
W0RK(DELTAI-1+I+(J-1)*N)  = DELTA(I,J)  = Au 
for  I = 1, . . . ,n,  and  J = 1, . . .,m. 

DELTAI  = 1 . 

> WORK(EPSI)  is  the  first  element  of  an  n x q array  EPS  containing  the  estimated 

errors  E in  the  response  variables  at  the  solution,  where 
W0RK(EPSI-1+I+(L-1)*N)  = EPSCI.L)  = 
for  I = 1, . . . , n,  and  L = 1, . . . , q. 

EPSI  = nm  + 1 . 

> WORK(XPLUSI)  is  the  first  element  of  an  n x m array  XPLUSD  containing  the  final 

estimates  of  the  explanatory  variable  X,  where 
W0RK(XPLUSI-1+I+(J-1)*N)  = XPLUSDCI,  J)  = Xu  = Xu  + Au 
for  I = 1, . . . , n,  and  J = 1, . . . , m. 

XPLUSI  = nm  + nq  + 1 . 
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> WORK(FNI)  is  the  iirst  element  of  an  n X 9 array  FN  containing  the  final  estimates 

of  the  response  variable  V,  where 

WORK(FNI-l+I+a-l)*N)  = FNCI.L)  = Fn,  = /ilC®!  + 
for  I = 1, . . . , n,  and  L = 1, . . . , g. 

FNI  = 2nm  + ng  + 1 . 

> WORK(SDI)  is  the  first  element  of  a p X 1 array  SD  containing  the  standard 

deviations  of  the  function  parameters  /?,  i.e.,  the  square  roots 
of  the  diagonal  entries  of  the  covariance  matrix,  where 

W0RK(SDI-1+K)  = SD(K)  = 

for  K = 1, . . . ,p.  The  standard  deviations  are  only  computed  when 
the  third  digit  of  JOB  is  less  than  or  equal  to  1.  (See  §2.6.ii,  subrou- 
tine argument  JOB,  and  §4.B.)  Rows  of  SD  corresponding  to  fixed 
elements  of  and  to  elements  dropped  because  they  induced  rank 
deficiency,  are  set  to  zero. 

SDI  = 2nm  -|-  2ng  -|- 1 . 


> WORK(VCVI)  is  the  first  element  of  a p x p array  VCV  containing  the  values  of 

the  covariance  matrix  of  the  parameters  prior  to  scaling  by  the 
residucd  variance,  where 
W0RK(VCVI-1+I+(J-1)*(NP))  = VCV(I.J)  = 
for  I = l,...,p  and  J = l,...,p.  The  covariance  matrix  is  only 
computed  when  the  third  digit  of  JOB  is  less  than  or  equal  to  1.  (See 
§2.Bii,  subroutine  argument  JOB,  and  §4.B.)  Rows  and  columns  of 
VCV  corresponding  to  fixed  elements  of  /?,  and  to  elements  dropped 
because  they  induced  rank  deficiency,  are  set  to  zero. 

VCVI  = 2nm  H-  2ng  -f-  p -f  1 . 

> WORK(RVARI)  is  the  estimated  residual  variance  = 30,6)1^1.  (See  §4.B.) 

RVARI  = 2nm  + 2nq  -f  p + p^  + 1 • 

> WORK(WSSI)  is  30,6).  (See  §1.A.) 

WSSI  = 2nm  -f-  2ng  -|-  p -|-  p^  -f  2 . 


> WORK(WSSDEI)  is  360,6).  (See  §1.A.) 

WSSDEI  = 2nm  -H  2ng  -1-  p + p^  -f  3 . 
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> WORK(WSSEPI) 

> WORK(RCONDI) 

> WORK(ETAI) 

WORK(OLMAVI) 

WORK(TAUI) 

WORK(ALPHAI) 

WORK(ACTRSI) 

WORK(PNROMI) 

WORK(RNORSI) 

WORK(PRERSI) 

WORK(PARTLI) 

WORK(SSTOLI) 

WORK(TAUFCI) 

WORK(EPSMAI) 

WORK(BETAOI) 


is  Se0,6).  (See  §1.A.) 

WSSEPI  = 2nm  + 2nq  + p + + 4 . 

is  K~^{MsJp)  at  the  solution.  (See  §4.C.) 

RCONDI  = 2nm  + 2nq  + p + p^  + 5 . 

is  the  relative  error  in  the  function  values  computed  within  FCN. 
ETAI  = 2nm  + 2nq  + p + p*  + 6 . 

is  the  average  number  of  steps  required  to  obtain  the  Levenberg- 
Marquardt  parameter. 

is  the  trust  region  radius  at  the  time  the  computations  stopped. 

is  the  Levenberg-Marquardt  parameter  at  the  time  the  computa- 
tions stopped. 

is  the  actual  relative  reduction  in  5(/3, 6)  from  the  last  iteration, 
is  the  norm  of  the  scaled  values  of  and  S. 
is 

is  the  predicted  relative  reduction  in  5(/3,^)  from  the  last  iteration. 

is  the  stopping  tolerance  used  to  detect  parameter  convergence. 

is  the  stopping  tolerance  used  to  detect  sum  of  squares  convergence. 

is  the  factor  used  to  compute  the  initial  trust  region  radius. 

is  machine  precision,  i.e.,  the  smallest  value  ^ such  that  1 -1-  ^ > 1 . 

is  the  first  element  of  a p x 1 array  BETAO  containing  the  initial 
estimates  of  the  function  parameters  /3°,  where 
WORK(BETAIO-1+K)  = BETAO (K)  = 
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WORK(BETACI) 


WORK(BETASI) 


WORK(BETANI) 


WORK(SI) 


WORK(SSI) 


WORK(SSFI) 


WORK(QRAUXI) 


for  K = 1, . . . ,p.  For  implicit  models,  BETAO  is  the  initial  value  used 
for  the  last  value  of  the  penalty  parameter. 

is  the  first  element  of  a p x 1 array  BETAC  containing  the  current 
working  estimates  of  the  p unfixed  function  parameters  where 
W0RK(BETACI-1+K)  = BETAC  (K)  = 
for  K = 1, . . . ,p. 

is  the  first  element  of  a p x 1 array  BETAS  containing  the  saved 
working  estimates  of  the  p unfixed  function  parameters  where 
W0RK(BETASI-1+K)  = BETAS (K)  = 
for  K = 1, . . . ,p. 

is  the  first  element  of  a px  1 axray  BET  AN  containing  the  new  working 
estimates  of  the  p unfixed  function  parameters  where 
W0RK(BETANI-1+K)  = BETAN(K)  = 
for  K = 1, . . . ,p. 

is  the  first  element  of  a p x 1 array  S containing  the  step  in  the  p 
unfixed  function  parameters  at  the  last  iteration,  where 
W0RK(SI-1+K)  = S(K) 
for  K = 1, . . . ,p. 

is  the  first  element  of  a p x 1 array  SS  containing  the  scale  of  the  p 
unfixed  function  parameters  where 
W0RK(SSI-1+K)  = SS(K)  = SCALE{4k} 
for  K = 1, . . . ,p. 

is  the  first  element  of  a p X 1 array  SSF  containing  the  scale  of  all 
of  the  function  parameters  /3,  where 
W0RK(SSFI-1+K)  = SSF(K)  = SCALE{/3k} 
for  K = 1, . . . ,p. 

is  the  first  element  of  a p X 1 array  QRAUX  used  during  the  compu- 
tations, where 

W0RK(QRAUXI-1+I)  = QRAUX(I) 
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for  I = 1, . . . ,p. 

WORK(U)  is  the  first  element  of  a p X 1 array  U used  during  the  computations, 
where 

W0RK(UI-1+I)  = U(I) 
for  I = 1, . . . ,p. 

WORK(FSI)  is  the  first  element  of  an  n X g array  FS  containing  the  saved  esti- 
mated errors  E*  m the  response  variable,  where 
W0RK(FSI-1+I+L*N)  = FS(I,L)  = 
for  I = 1, . . . , n,  and  L = 1, . . . , 9. 


WORK(FJACBI)  is  the  first  element  of  an  nxpxq  array  FJACB  containing  the  weighted 
partial  derivative  with  respect  to  the  p unfixed  function  parameters 
P,  where 

WORK  (FJACBI- 1+1+ (K-l)*N+a-l)*N*NP) 

= FJACB(I,K,L)  = 

for  I = 1,  ...,n,  K = 1,  ...,p,  L = and  Wi'^Wi  = rtfej. 

The  derivatives  are  the  values  evaluated  at  the  beginning  of  the  last 
iteration  unless  the  user  requested  that  the  covariance  matrix  be 
computed  using  the  final  solution,  in  which  case  they  are  the  values 
obtained  at  the  final  solution.  (See  §2.BJi,  subroutine  argument 
JOB.) 

WORK(WEII)  is  the  first  element  of  an  LDWE  X LD2WE  X g array  WEI  containing  the 
Cholesky  factorizations  for  the  weights  Wc  specified  in  WE,  where 
W0RK(WE1I-1+I+(L1-1)*LDWE+(L2-1)*LDWE*LD2WE) 

= WE1I(I,L1,L2) 

for  I = 1,...,LDWE,  LI  = 1,...,LD2WE,  and  L2  = l,...,g.  WEII 
specifies  the  factorization  in  the  same  manner  that  subroutine  ar- 
gument WE  is  used  to  specify  We.  (See  §2.B.ii.) 

WORK(DIFFI)  is  the  first  element  of  an  g X (p  -|-  m)  array  DIFF  containing  the 
relative  differences  between  the  user  supplied  derivatives  and  the 
finite  difference  values  they  were  checked  against,  where 
W0RK(DIFFI-1+L+(J-1)+NQ)  = DIFFCL.J) 
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for  L = 1, . . . , g,  and  J = 1, . . . , (p  + m). 

WORK(DELTSI)  is  the  first  element  of  an  n x m array  DELTAS  containing  the  saved 
working  estimates  of  the  errors  A'  in  the  explanatory  variables, 
where 

WORK  (DELTSI- 1+1+ (J-1)*N)  = DELTASCl,  J)  = Aj;, 

for  I = 1, . . . , n,  and  J = 1, . . . , m.  If  ISODR  is  false,  then  this  array 

is  equivalenced  to  the  array  DELTA  beginning  in  WORK(DELTAI). 

WORK(DELTNI)  is  the  first  element  of  an  n x m array  DELTAN  containing  the  new 
working  estimates  of  the  errors  A”  in  the  explanatory  variables, 
where 

W0RK(DELTNI-1+I+(J-1)*N)  = DELTAN(I,J)  = AJJ 

for  I = 1, . . . , n,  and  J = 1, . . . , m.  If  ISODR  is  false,  then  this  array 

is  equivalenced  to  the  array  DELTA  beginning  in  WORK(DELTAI). 


WORK(TI)  is  the  first  element  of  am  n X m array  T used  in  the  computations, 
where 

WORK (TI- 1+1+ (J-1)*N)  = TCI.J) 

for  I = 1, . . . , n,  and  J = 1, , . . , m.  If  ISODR  is  false,  then  this  array 
is  equivalenced  to  the  array  DELTA  beginning  in  WORK(DELTAI). 

WORK(TTI)  is  the  first  element  of  an  n x m array  TT  containing  the  scale  of  each 
the  estimated  errors  A in  the  explanatory  variable,  where 
W0RK(TTI-1+I+(J-1)*N)  = TT(I,J)  = scale{Au} 
for  I = 1, . . . , n,  and  J = 1, . . . , m.  If  ISODR  is  false,  then  this  array 
is  equivalenced  to  the  array  DELTA  beginning  in  WORK(DELTAI). 

WORK(OMEGAI)  is  the  first  element  of  a q x q array  OMEGA  used  during  the  compu- 
tations, where 

W0RK(0MEGAI-l+Ll+a2-l)*NQ)  = 0MEGAai,L2) 

for  LI  = 1, . . . , q,  and  L2  = 1, . . . , q.  If  ISODR  is  false,  then  this  array 

is  equivalenced  to  the  array  DELTA  beginning  in  WORK(DELTAI). 


WORK(FJACDI)  is  the  first  element  of  an  n X m X g array  FJACD  containing  the 
weighted  partial  derivative  with  respect  to  A,  where 
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WORK(FJACDI-l+I+(J-l)*N+a-l)*N*M) 

= FJACD(I,J,L)  = 

for  I = J = L = and  Wi'^Wi  = uti- 

The  derivatives  are  the  values  evaluated  at  the  beginning  of  the  last 
iteration  unless  the  user  requested  that  the  covariance  matrix  be 
computed  using  the  final  solution,  in  which  case  they  are  the  values 
obtained  at  the  final  solution.  (See  §2.Bii,  subroutine  argument 
JOB.)  If  ISODR  is  false,  then  this  array  is  eqidvalenced  to  the  array 
DELTA  beginning  in  WORK(DELTAI). 

WORK(WRKII)  is  the  first  element  of  an  n X m X g array  WRKl  required  for  work 
space,  where 

WORK(WRKlI-l+I+(J-l)*N+a-l)*N*NQ)  = WRK1(I,J,L) 

for  I = 1, . . . , n,  J = 1, . . . , m,  and  L = 1, . . . , g.  If  ISODR  is  false, 

then  this  array  is  equivalenced  to  the  array  DELTA  beginning  in 

WORK(DELTAI). 

W0RK(WRK2I)  is  the  first  element  of  an  n x g array  WRK2  required  for  work  space, 
where 

W0RK(WRK2I-1+I+(L-1)*N)  = WRK2(I,L) 
for  I = 1, . ..,71,  and  L = 1,.  ..,g. 

W0RK(WRK3)  is  the  first  element  of  a p x 1 array  WRK3  required  for  work  space, 
where 

W0RK(WRK3-1+K)  = WRK3(K) 
for  K = 1, . . . ,p. 

W0RK(WRK4I)  is  the  first  element  of  an  m x m array  WRK4  required  for  work  space, 
where 

W0RK(WRK4I-1+J1+(J2-1)*M)  = WRK4(J1,J2) 
for  J1  = 1, ...  ,771,  and  J2  = 1, . . .,7n. 

W0RK(WRK5I)  is  the  first  element  of  an  tti  array  WRK5  required  for  work  space, 
where 

W0RK(WRK5I-1+J)  = WRK5(J) 
for  J = 1,  ...  ,771. 
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W0RK(WRK6I)  is  the  first  element  of  an  n x p X g array  WRK6  required  for  work 
space,  where 

W0RK(WRK6I-1+I+(K-1)*N+(L-1)*N*NQ)  = WRK6(I,K,L) 
for  I = 1,  K = and  L = 

W0RK(WRK7I)  is  the  first  element  of  an  5g  array  WRK7  required  for  work  space, 
where 

W0RK(WIIK7I-1+J)  = WRK7(J) 
for  J = 


5.B.  Extracting  Information  from  Vector  IWORK 

Upon  return  from  a call  to  ODRPACK,  array  IWORK  contains  values  that  may  be  of 
interest  to  the  user.  To  extract  information  from  IWORK,  the  following  declaration 
statement  must  be  added  to  the  user’s  program 

INTEGER 

+ MSGBI , MSGDI , IFIX2I , ISTOPI , 

+ NNZWI.NPPI.IDFI, 

+ JOBI,IPRINI,LUNERI,LUNRPI, 

+ NROWI.NTOLI.NETAI, 

+ MAXITI , NITERI . NFEVI , N JEVI , INT2I , IRANKI , LDTTI , 

+ LIWKMN 

where  MSGBI  through  LDTTI  are  variables  that  indicate  the  starting  locations  within 
IWORK  of  the  stored  values,  and  LIWKMN  is  the  minimum  acceptable  length  of  array  IWORK. 
The  appropriate  values  of  MSGBI  through  LDTTI  are  obtained  by  invoking  subroutine 
SIWINF  when  using  either  of  the  single  precision  ODRPACK  subroutines  SODR  or  SODRC, 
and  by  invoking  DIWINF  when  using  either  of  the  double  precision  subroutines  DODR  or 
DODRC.  The  call  statements  for  SIWINF  and  DIWINF  have  the  Scune  argument  lists.  To 
invoke  either  subroutine,  use 

CALL  <iwinf> 

+ (M.NP.NQ, 

+ MSGBI, MSGDI, IFIX2I, ISTOPI. 

+ NNZWI.NPPI.IDFI, 

+ JOBI.IPRINI.LUNERI.LUNRPI, 

+ NROWI.NTOLI.NETAI, 
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+ MAXITI , NITERI , NFEVI , N JEVI , INT2I , IRANKI , LDTTI , 

+ LIHKMN) 

where  SIWINF  should  be  substituted  for  <ivinf>  when  using  SODR  and  SODRC,  and 
DIWINF  should  be  substituted  for  <iwinf  > when  using  DODR  and  DODRC.  Note  that  the 
values  of  N,  NP,  and  NQ  must  be  input  to  SIWINF  and  DIWINF  with  exactly  the  same 
values  as  were  used  in  the  original  call  to  ODRPACK. 

In  the  following  descriptions  of  the  information  returned  in  IWORK,  > indicates  values 
that  are  most  likely  to  be  of  interest. 

> IWORK (MSGBI)  is  the  first  element  of  a 1 + (g  X p)  array  MSGB  used  to  indicate 
the  results  of  checking  the  partial  derivatives  with  respect  to  ^ at 
observation  NROW.  (See  IWORK(NROWI)  below.) 

The  value  of  IWORK  (MSGBI)  summarizes  the  results  over  all  0. 

• If  IWORK (MSGBI)  < 0 then 

the  partial  derivatives  with  respect  to  were  not  checked. 

• If  IWORK (MSGBI)  = 0 then 

the  partial  derivatives  with  respect  to  each  of  the  ;9k>  K = 1, . . . , p, 
for  each  of  the  q responses  appear  to  be  correct. 

• If  IWORK (MSGBI)  = 1 then 

the  partial  derivative  with  respect  to  at  least  one  of  the  /?k,  K = 

1. .  . .,p,  appears  to  be  questionable  for  at  least  one  of  the  q re- 
sponses. 

• K IWORK (MSGBI)  = 2 then 

the  partial  derivative  with  respect  to  at  least  one  of  the  /3k.  K = 

1..  ..,p,  appears  to  be  seriously  questionable  for  at  least  one  of 
the  q responses. 

The  value ofIWORK(MSGBI+L+(K-l)*NQ),  L = 1,...,9,  K = l,...,p, 
indicates  the  individual  results  of  checking  the  partial  derivative  of 
the  Lth  response  with  respect  to  each  /3k,  where  for  L = 
and  K = 1, . . . ,p  : 

• If  IWORK  (MSGBI+L+(K-1)*NQ)  = -1  then 

the  partial  derivative  of  the  Lth  response  with  respect  to  (3^.  was 
not  checked  because  (3k  was  fixed. 

• If  IWORK (MSGBI+L+(K-1)*NQ)  = 0 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  /3k  ap- 
pears to  be  correct,  i.e.,  the  relative  difference  between  its  value 
and  the  finite  difference  approximation  it  is  checked  against  is 
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within  the  required  tolerance. 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 1 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  is 
questionable  because  the  user  supplied  derivative  and  the  finite 
difference  value  it  is  checked  against  are  both  zero. 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 2 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  ydx  is 
questionable  because  either  the  user  supplied  derivative  is  ex- 
actly zero  and  the  finite  difference  value  it  is  checked  against  is 
only  approximately  zero,  or  the  user  supplied  derivative  only  ap- 
proximately zero  and  and  the  finite  difference  value  it  is  checked 
against  is  exactly  zero. 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 3 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  /3k  is 
questionable  because  either  the  user  supplied  derivative  is  exactly 
zero  and  the  finite  difference  value  it  is  checked  against  is  not 
even  approximately  zero,  or  the  user  supplied  derivative  not  even 
approximately  zero  and  and  the  finite  difference  value  it  is  checked 
against  is  exactly  zero. 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 4 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  /3k  is 
questionable  because  the  finite  difference  value  it  is  being  checked 
against  is  questionable  due  to  a high  ratio  of  relative  curvature 
to  relative  slope  or  to  an  incorrect  scale  value. 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 5 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  /3k  is 
questionable  because  the  finite  difference  value  it  is  being  checked 
against  is  questionable  due  to  a high  ratio  of  relative  curvature 
to  relative  slope. 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 6 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  /3k 
is  questionable  because  it  does  not  agree  with  the  finite  dif- 
ference value  it  is  being  checked  against  to  the  required  toler- 
ance, although  the  values  do  agree  in  their  first  two  digits.  (See 
IWORK(NTOLI)  below.) 

• If  IW0RK(MSGBI+L+(K-1)*NQ)  = 7 then 

the  partial  derivative  of  the  Lth  response  with  respect  to  /3k  is 
seriously  questionable  because  it  has  fewer  than  two  digits  agree- 
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ment  with  the  finite  difference  value  it  is  being  checked  against. 
MSGBI  = 1 . 

> IWORK(MSGDI)  is  the  first  element  of  a 1 + (g  X m)  array  MSGD  used  to  indicate 

the  results  of  checking  the  partial  derivatives  with  respect  to  A, 
analagous  to  MSGB  described  above.  The  values  in  MSGD  have  the 
same  meanings  as  those  used  to  indicate  the  results  of  checking 
the  partial  derivatives  with  respect  to  /3,  except  that  the  value  of 
IWORK(NSGDI)  summarizes  the  results  over  all  columns  of  A,  and  the 
values  of  IW0RK(MSGDI+L+(J-1)*NQ),  L = J = l,...,m, 

indicates  the  individual  results  for  checking  the  partial  derivative  of 
the  Lth  response  with  respect  to  the  Jth  column  of  A at  observation 
NROW.  (See  IWORK(NROWI)  below.) 

MSGDI  = qp  + 2 . 

> IW0RK(IFIX2I)  is  the  first  element  of  a p X 1 array  IFIX2  containing  values  used 

to  indicate  whether  a given  parameter  is  unfixed,  fixed  or  dropped 
because  it  induced  rank  deficiency,  where 
W0RK(IFIX2I-1+K)  = IFIX2(K) 
for  K = 1, . . . ,p.  If 

• IFIX2(K)  = 1 the  parameter  was  unfixed, 

• IFIX2(K)  = 0 the  parameter  was  fixed, 

• IFIX2(K)  = — 1 the  parameter  was  dropped,  and 

• IFIX2(K)  = —2  no  parameters  were  estimated. 

IFIX2I  = qp  + qm  + 3 . 

> IWORK(ISTOPI)  is  value  of  ISTOP  returned  from  the  last  call  to  subroutine  FCN. 

ISTOPI  = qp  + qm  + p + 3 . 

IWORK(NNZWI)  is  the  number  of  nonzero  c error  weights,  tWej,  I = 1, . . . , n. 

IWORK(NPPI)  is  the  number  p of  function  parameters  actually  being  estimated, 
i.e.,  the  number  of  unfixed  parameters. 


IWORK(IDFI)  is  the  degrees  of  freedom,  p,  of  the  fit,  i.e.,  the  number  of  observa- 
tions with  nonzero  weights  minus  the  number  of  parameters  actually 
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IWORK(JOBI) 

IWORKCIPRINI) 

IWORK(LUNERI) 

IWORK(LUNRPI) 

IWORK(NROWI) 

IWORK(NTOLI) 

IWORK(NETAI) 

IWORK(MAXITI) 

IWORK(NITERI) 

IWORK(NFEVI) 

IWORK(NJEVI) 

IW0RK(INT2I) 

IWORK(IRANKI) 

IWORK(LDTTI) 


being  estimated. 

is  the  value  used  to  specify  problem  initialization  and  computational 
methods. 

is  the  print  control  value  used. 

is  the  logical  unit  number  used  for  error  reports. 

is  the  logical  unit  number  used  for  computation  reports. 

is  the  observation  NROW  at  which  the  derivatives  were  checked. 

is  the  number  of  digits  of  agreement  required  between  the  numerical 
derivatives  and  the  user  supplied  derivatives  for  the  user  supplied 
derivatives  to  be  considered  correct. 

is  the  number  of  good  digits  in  the  function  results  returned  by  user 
supplied  subroutine  FCN. 

is  the  maximum  number  of  iterations  allowed. 

is  the  number  of  iterations  taken. 

is  the  number  of  function  evaluations  made. 

is  the  number  of  Jacobian  matrix  evaluations  made. 

is  the  number  of  internal  doubling  steps  taken  at  the  time  the  com- 
putations stopped. 

is  the  rank  deficiency  at  the  solution. 

is  the  leading  dimension  of  the  n x m array  TT.  (See  §5.A.) 
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